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I .  Affine  Transformations  and  Tracking . 

When  a  dynamic  three-dimensional  scene  is  observed  via  an 
optical  projection,  a  quite  complex  class  of  motions  are  induced 
in  the  image  plane  [5, 6, 7, 9].  In  general,  this  class  of  motions 
is  highly  non-linear,  being  dependent  on  the  geometry  of  the 
objects  being  observed  as  well  as  their  trajectories  in  space [3]. 
nevertheless ,  in  many  cases  the  motion  is  approximated  closely 
by  translation,  magnification  and  rotation  in  the  image  plane. 

It  is  easy  to  see  that  this  approximation  is  best  for  motion  in 
space  which  consists  of  translation  and  rotation  about  a  line 
parallel  to  the  bare  sight. 

A  better  approximation  results  by  consideration  of  the  full 
affine  group  in  the  plane,  which  includes  shearing  in  two  direct¬ 
ions  as  well  as  the  motions  mentioned  above.  By  definition,  an 
affine  transformation  in  the  plane  R  is  of  the  form 

T(y)  =  Ay  +  a,  yeR2  (1.1) 

2 

where  A  is  a  non-singular  2x2  matrix  and  asR  is  considered  as  a 
column  vector  [4,10].  The  set  of  all  such  transformations  T  is 
called  the  general  affine  group  and  is  denoted  GA(2) .  It  is 
easily  seen  that  the  subset  consisting  of  translations,  magni¬ 
fications  and  rotations  forms  a  subgroup,  which  we  denote  by 
SA(2).  In  order  that  T(y)  =  Ay  +  a  belong  to  SA(2)  it  is  nec¬ 
essary  and  sufficient  that  A^  =  an^  A^  =  -A21‘  In  th-'-s 

2  2  1/2 

case,  the  magnification  factor  is  (A^  +  A21^  and  t*ie  rotati°n 
angle  is  atn (A^/A^) . 

In  order  to  consider  dynamic  images,  it  is  necessary  to  allow 


2 


A  and  a  in  (1.1)  to  depend  on  time.  This  gives  rise  to  a  tra- 

2 

jectory  u(t,y)  for  each  yeR  given  by 

u(t,y)  =  A ( t) y  +  a ( t)  (1.2) 

where  (A(t) ,  a(t))  e  GA(2),  and  we  require  that  A(0)=I,  a(0)=  0 
in  order  that  the  trajectory  pass  through  y  at  tine  t=0;  i.e., 
u (0 ,y)  =  y . 

As  in  [4],  though  only  for  linear  transformations,  we  may 
realize  the  pair  (A(t),  a(t))  as  the  solution  of  a  linear  system 
of  differential  equations.  Let  us  define 

A  (t)  =  A  ( t)  A~ 1  ( t )  (1.3a) 


A  ( t)  =  a  (t)  -  A  (t)  a  (t)  ,  (1.3b) 

from  which, 

A  ( t)  =  A  ( t)  A  ( t)  ,  A  ( 0)  =  I  (1.4a) 


a (t)  *  A ( t)  +  A ( t)  a ( t) ,  a(0)  =  0  (1.4b) 

We  may  summarize  the  correspondences  defined  by  (1.3)  and  (1.4) 
as  follows: 


Theorem  1.1:  Equations  (1.3)  and  (1.4)  establish  a  one-to-one 
correspondence  between  differentiable  curves  (A(t),  a(t))  in 
GA (2)  satisfying  A(0)  =  I,  a(0)  =  0  and  continuous  curves 

2 

(A(t),  \(t))  where  A (t)  is  an  arbitrary  2x2  matrix  and  A (t)  eR  . 
Moreover,  in  order  that  (A, a)  belong  to  SA(2)  it  is  necessary 


and  sufficient  that  A^  =  A22  and  A^  =  ~^21' 


The  first  part  of  the  above  theorem  apparent  from  (1.3)  and 
(1.4).  A  rigorous  and  detailed  proof  proceeds  exactly  as  given 
in  [4]  for  linear  transformations.  The  last  part  can  be  deduced 
by  a  few  calculations  using  the  fact  that  elements  of  SA(2)  satisfy 

A11  "  A22  and  A12  *  -A21' 


3 


Now,  if  we  differentiate  (1.2)  with  respect  to  t,  use  (1.4) 
and  (1.2)  again,  we  obtain 

— |^— (t,y)  =  A ( t)  u(t,y)  +  \(t).  (1.5) 

Of  course,  v(t,y)  =  —  (t,y)  is  the  velocity  field  along  the 

J  t 


trajectory  u(t,y).  Equ.  (1.5)  shows  that  the  velocity  at  a  point 
u  depends  on  u  as  well  as  t  and  is  therefore  not  spatially  in¬ 
variant  . 

Note  that  (1.5)  in  fact  gives  the  differential  equation  for 
an  arbitrary  affine  trajectory,  and  when  A(t)  is  restricted  as 
in  Theorem  1.1,  it  gives  the  equation  for  a  trajectory  under  the 
restricted  group  of  motions  SA(2) .  By  virtue  of  (1.2) ,  we  see 
that  ( A ( t ) ,  a(t))  obtained  from  (1.4)  may  be  considered  as  a 
fundamental  system  of  solutions  to  the  evolution  equation  (1.5). 
Now  it  is  important  to  note  that  the  fundamental  system  of  solu¬ 
tions  is  completely  determined  by  the  pair  (A(t),  A  ( t ) )  which 
is  spatially  invariant,  being  a  function  of  time  only.  To  es¬ 
tablish  convenient  terminology,  let  us  give  the  following 
Definition  1.1:  The  pair  { A ( t ) ,  A(t))  is  the  generalized  velocity 
field  of  the  family  of  affine  trajectories  u(t,y)  defined  by  (1.5). 


We  may  now  state 

Theorem  1.2:  A  family  u(t,y)  of  affine  trajectories  satisfying 
u(0,y)  =  y  is  completely  determined  from  its  generalized  velocity 
field,  which  is  spatially  invariant,  via  (1.4)  and  (1.2).  More¬ 
over,  the  absolute  velocity  v  at  a  point  u  on  a  trajectory  is 
given,  as  in  (1.5),  by  v  =  A (t) u  +  X  (t)  . 


rh 

L:'2j 


I 


Let  us  write  A(t) 


,  A  (t) 


and  expand  the  equa- 
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tion  u  =  A  u  +  X  (where  the  t-dependence  has  been  suppressed)  in 
the  form 
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In  this  way  we  can  identify  individual  vector  fields  v^(u),..., 
v^ (u)  and  write  (1.6)  in  the  form 

-g§-  =  l  ^(t)  V1(u)  (1.7) 

i=l 


In  a  similar  manner  for  SA(2) ,  we  write 


so  that 
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(1.85 


1  4 

allowing  four  vector  fields  v  (u),...,  v  (u)  to  be  inaentified. 
Finally,  we  rewrite  (1.8)  as 


3u 


4 

l  X  (t)  v1  (u) 
i=l 


(1.9) 


It  should  be  noted  that  the  functions  v1  defined  by  (1.7)  or 
(1.9)  are  characteristic  of  the  class  of  motions  under  consider¬ 
ation,  and  more  general  classes  of  motion  can  be  treated  by 

1  2 

consideration  of  other  generators  v  ,v  , '  '  ' 


I 


.  In  the  cases  of 


5 


interest,  the  sets  of  vector  fields  derived  above  define  the  Lie 
algebras  [2]  of  GA(2)  and  SA(2)  .  Any  vector  field  v:  R4-*-R“  induces  a 
differential  operator  Y  ,  called  an  infinitesimal  transformation, 
which  is  defined  by 


Y 

v 


V, 


1 


+ 


V, 


(1.10) 


where  v^(y)  and  v^(y)  are  the  components  of  v(y).  In  Tables 
1  and  2  we  list  the  infinitesimal  transformations  for  the  groups 
GA ( 2 )  and  SA(2) ,  given  in  terms  of  a  variable  x=(x1,x2)  for  later 
application 

Table  1.  Infinitesimal  transformations  for  GA(2) . 

X1  =  Jx^  X3  =  X1  Jx^  X5  =  X2  3x^ 

X2  =  X4  =  X1  JxJ  X6  =  X2 


Table  2. 


X 


1 


X 


2 


Infinitesimal  transformations  for  SA(2) 


3 

Y 

3  3 

3X1 

X3 

X1  axL  x2  3x2 

3 

X4 

3  3 

3x2 

X1  3x2  X2  3x^ 

Note  that  u(t,y)  given  by  (1.2)  may  be  regarded  as  the  lo¬ 
cation  at  time  t  of  the  particle  which  was  at  y  at  time  t=0. 

An  observar  at  some  point  x  will  observe  this  particle  provided 
that  x  =  u(t,y)  =  A(t)y  +  A ( t ) .  We  may  solve  this  equation  for 

y  to  obtain  y  =  A  ^(t)  (x  -  a(t))  .  Thus  we  define  the  trace  of 
2 

the  point  x  R  to  be 


-1  2 
s(t,x)  =  A  (t)(x-a(t)),  xeR  . 
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(1.11) 


We  may  interpret  s(t,x)  to  be  the  particle  which  will  arrive  at 
x  at  time  t. 

Let  us  now  consider  a  two-dimensional  image,  represented  by 
2 

a  runction  r:R  -►R,  and  suppose  that  the  image  r  is  subjected  to 
an  affine  transformation  (A(t),  a(t)).  Here  a  value  f(y)  is  re¬ 
garded  as  a  feature  which  propogates  along  the  trajectories  of 
the  motion.  This  is  an  extremely  powerful,  and  somewhat  re¬ 
strictive,  assumption  which  is  not  always  valid  in  real  images. 

For  example,  it  is  violated  by  changes  in  radiance  values  which 
vary  as  a  function  of  the  angle  of  incident  illumination.  On 
the  other  hand,  it  is  valid  in  most  instances  over  short  time 
intervals  and  deviations  from  this  assumption  may  frequently  be 
treated  as  higher  order  effects. 

In  any  event,  if  the  feature  f (y)  is  propagated  along  tra¬ 
jectories,  then  a  stationary  observor,  say  at  point  x,  will  ob¬ 
serve  a  value  F(t,x)  =  f(s(t,x))  at  time  t,  since  s(t,x)  re¬ 
presents  the  particle  arriving  at  x  at  time  t.  We  may  now  state 
a  most  important  result. 

Theorem  1.3:  Let  a  time- varying  image  F  be  given  by 

F ( t , x)  =  f (s ( t , x ) )  (1.12) 

where  s(t,x)  is  an  affine  trace  as  in  (1.11)  with  generalized 
velocities 


7 

-F  5 

- 4  =  >  \.  (t)  X.  F  ,  (1.13) 

i=l  1  1 

where  X, are  given  in  Table  1. 

1  O 

A  similar  result  holds  if  we  restrict  to  SA(2),  using 
X^,...,X^  from  Table  2. 


Proof:  This  result  may  be  deduced  from  results  given  in 

[7],  provided  we  compensate  for  the  change  from  "left  invariance" 
in  that  development  to  the  "right  invariance"  of  the  current 
treatment.  However,  a  direct  proof  is  instructional  and  will 
be  outlined  herein.  We  first  show 
Lemma  1.3.1:  For  s(t,x)  given  by  (1.11) 


IF"  -Jx  Ai(t)xis- 


(1.14) 


Proof:  First  note  that  by  direct  calculation  we  have  E'4Xis= 

V  \  .  X  .  A- 1  ( x-  a )  =  A_L(y.\.X.x)  =  A-1  ( )  A  .v1  (x)  )  =A_1(Ax+X).  Also, 

‘-ii  -li  ‘-i 

noting  that  — ~  A  ^  =  -A  ^AA  ^ ,  we  have  — ~  A  (x~a)  = 

-A~IAA-1 (x-a) -A_1a  =  -A_1A(x-a)  -  A_L(\  +  Aa)  =  -A~1(Ax  +  \). 

Hence,  the  desired  result  follows. 

Returning  to  the  proof  of  Theorem  1.3,  we  have 

y  X.X.F(tx)  =  l  X.(t)vFx)  l—f  (s(t,x)  ) 

1  1  -i  1  3 
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=  y  x . (t)  x . s_  ( t ,x) 

T  -  IK 


3f 

3sk  (s) 


3s 
3 1 


k  ,,  .  3 f 

sk 


Is) 


3f (s (t,x) ) 
2t 


3F  ,  ... 

=  -  * —  (t,x)  ,  as  desirec. 

<3  t 

Theorem  1.3  appears  to  be  fundamental  to  the  analysis  of 
motion  in  dynamic  images.  As  is  evident  from  the  proof,  an 
analogue  is  valid  in  a  much  more  general  setting.  In  fact, 
scrutiny  of  the  proof  shows  that  it  depends  mainly  on  Lemma  1.3.1. 
Consequently,  the  theorem  will  hold  for  any  class  of  motions  for 
which  a  suitable  form  of  the  "trace”  lemma  can  be  obtained.  The 
■ignj f icance  of  Theorem  1.3  lies  chiefly  in  the  fact  that  the 
generalized  velocities  (which  are  usually  unknown)  appear  as 
linear  coefficients  in  (1.13),  along  with  quantities  which  can 
be  calculated  from  the  data  F(t,x) . 

The  main  problem  with  the  extraction  of  the  generalized 
velocities  from  (1.13)  is  the  general  lack  of  numerical  pre¬ 
cision  in  the  calculation  of  the  derivatives  from  real  data 
(e.g.,  digitized  video).  In  subsequent  sections  we  shall  show 
how  to  incorporate  (1.13)  in  a  feedback  loop  which  is  very  stable 
and  how  to  obtain  an  equivalent  formulation  based  on  integration 


rather  than  differentiation. 
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II.  A  Velocity  Feedback  Tracker 

The  theoretical  results  of  this  section  result  from  research 

* 

done  under  a  separate  contract  and  for  which  a  publication  is 
in  preparation.  In  view  of  the  fact  that  the  techniques  have 
been  incorporated  in  the  experimental  portion  of  this  report, 
these  results  will  be  presented  in  this  report  in  the  context 
of  affine  transformations. 

Let  absolute  image  coordinates  be  denoted  by  y  =  (y1»y2) 
and  introduce  additional  coordinates  as  follow:  Let  coordinates 
z  be  established  relative  to  a  moving  target,  and  let  coordinates 
x  be  established  in  a  movable  "window".  We  assume  that  the 
motions  of  both  the  target  and  the  window  may  be  described  by 
affine  transformations  relative  to  absolute  image  coordinates. 

It  is  assumed  that  the  motion  of  the  window  may  be  chosen  at 
will,  while  the  motion  of  the  target  is  prescribed  (e.g.  by 
nature)  and  is  unknown. 

3y  the  affine  assumption,  we  may  describe  the  transformation 
from  window  coordinates  to  image  coordinates  by 

yw(t,x)  =  A ( t) x  +  a(t)  ,  xeR*  ,  (2.1) 

where  (A(t) ,  a(t))  is  a  suitable  family  of  affine  transformations. 
Similarly,  the  transformation  from  target  coordinates  to  image 
coordinates  is 

r2 

yT(t,z)  =  B  ( t)  z  +b(t)  ,  ze  T  (2.2) 

for  suitable  (B(t),  b(t))  e  GA{2).  Let  us  denote  the  respective 
generalized  velocity  fields  by  (Aft,  X^)  and  (A \Q) . 

By  equating  y^,(t,z)  =  yw(t,x)  we  may  solve  for  the  point  z 
*  OMR  Contract  MOO 14-76-C-1136 
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on  the  target  which  arrives  at  point  x  in  the  window  at  time  t, 
to  obtain: 


z(t,x)  =  3~ 1 (Ax  +  a  -  b) ,  (2.3) 

where  dependence  on  t  has  been  suppressed  on  the  right.  We  note 
that  z(t,x)  in  (2.3)  may  be  regarded  as  a  trace  in  the  sense 
of  the  previous  section.  By  an  application  of  Lemma  2.3.1  we 
have : 

Theorem  2.1:  There  exists  a  generalized  velocity  field  (F(t),y(t)) 
such  that 


3z  ( t ,  x ) 
3t 


6 

■'£"  Y-  (t)  X .  z  (t  ,x)  , 
i=l  1  1 


(2.4) 


where  y(t) 


1 

-< 

CJ 

ui 

_ 1 

i 

-< 

I 

r  1  = 

_Y4  Y6_ 

and  the  operators  X^,...,Xg  are 


given  in  window  coordinates  as  in  Table  1. 


Now,  let  f(z)  be  a  feature  of  the  target,  measured  at  point 
z,  and  assume  that  this  feature  propogates  along  the  target 
tro jectories .  An  observor  at  point  x  in  the  window  therefore 
observes  data  F(t,x)  =  f(z(t,x)),  inasmuch  as  z(t,x)  is  the  point 
which  arrives  at  x  at  time  t.  From  Theorem  1.3  we  obtain 
Theorem  2.2:  In  the  above  context, 

6 

-  ~  (t,x)  =  y  Y  ■  ( t)  X  .  F  ( t ,  x)  (2.5) 

jt  i=l  1  1 


In  principle,  (2.5)  allows  the  determination  of  the  general¬ 
ized  velocities  of  the  target  relative  to  the  window.  Since 
we  have  free  choice  of  window  velocities,  relative  to  the  image 
coordinates,  this  is  tantamount  to  measurement  of  absolute  target 


11 


velocities.  The  conversion  process  will  now  be  described. 

Let  us  denote  by  X  the  window  space,  Y  the  image  space,  and 

let  T ( X)  and  T  (Y)  be  the  respective  tangent  spaces.  Since  the  map 

from  X  to  Y  is  y  =  A  x  +  a,  as  in  (2.1) ,  it  follows  that  the  in- 

★  ★ 

duced  map  on  the  tangent  spaces  is  simply  y  =  A  x  [1,2]. 

Now  the  velocity  field  ( T ,  y)  defines  a  map  X  -*■  T(X)  given  by 
* 

x  =  F  x  +  y  (see  (1.5)).  Accordingly,  a  velocity  field  (7,  X) 
is  induced  on  y.  which  maps  Y  -*■  T  (Y)  in  a  similar  fashion.  The 
velocity  field  (A,  X)  is  defined  by  the  commutative  diagram 


* 

Thus,  we  calculate  y  =  A  y  +  X ,  by  inverting  (A, a)  and  taking 

*  -1  -1 

the  upper  path,  to  be  given  by  y  =  A(F  A  (y-a)+y)  =  ArA  y  + 
Ay-ArA  ^a.  By  comparison,  we  obtain 

A  =  ArA-1  (2.6a) 

X  =  Ay-AFA^a  (2.6b) 

Now,  since  the  velocity  field  (F,y)  represents  the  difference  be¬ 
tween  target  and  window  velocities  in  the  window  coordinate  system, 
we  see  that  ( A , X )  must  represent  this  same  difference  relative  to 
the  absolute  image  coordinate  system.  That  is, 


A  -  ab  “aa 


(2.7a) 


X 


(2.7b) 
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We  may  summarize  these  results  in  a  useable  formas  follows: 
Theorem  2.3:  Let  (2.1)  and  (2.2)  define  the  motion  of  a  window 
and  a  target,  respectively,  relative  to  a  system  of  absolute  image 
coordinates,  and  let  (AA,  \^)  and  (A  ,  AQ)  be  the  corresponding 
generalized  velocities.  Further,  let  (F,y)  be  the  generalized 
velocities  of  the  target  relative  to  the  window,  as  determined  by 
(2.5)  . 

Then 

Ab  -  aa  =  AFA-1  (2.8a) 

A  -  Aa  =  Ay-AFA^a.  (2.8b) 

The  previous  theorem  immediately  suggests  an  algorithm  for 
determination  of  velocities  in  an  image.  More  generally,  the 
algorithm  performs  tracking  since,  as  will  be  seen,  the  result  is 
to  force  the  window  to  follow  the  target  by  emulation  of  velocities. 
The  algorithm  is  as  follows: 

Step  1.  Initialize  the  window  by  choice  of  A(0),  a(0).  In  the 
absence  of  a  priori  information,  initialize  AA(0)=0,  Aa(0)=0. 

Sample  window  values  F(tQ,x)  at  time  tQ=0. 

Step  2.  Sample  window  values  F(t  ,x)  at  time  t  =t  .+5.  Appro- 
^ n  n  n-i 

ximate  and  X^  F  at  various  points  in  the  window  and  form  a 
system  of  linear  equations  using  (2.5). 

Step  3.  Solve  the  resulting  linear  equations  for 

Step  4.  Replace  AA“-AA  +  AFA  ^  and  A  «-AA  +  Ay  -  AFA  1a. 

Note:  If  the  calculation  of  y.  y,....  were  exact,  this 

1.  r  Z 

would  result  in  A,-<-A_(t)  and  \.*-A-,(t  ) 

A  8  n  Aon 

Step  5.  Take  a  5  time  step  in  the  numerical  solution  A  =  AAA, 

a  =  y.  +  A, a  to  obtain  (Aft  ),  a  (t  )).  This  effectively  moves 
a  A  n  n 


the  window. 


Step  6.  Repeat  from  step  2. 

Emperical  results  indicate  that  the  above  algorithm  conveys 
rapidly  over  a  fairly  broad  range  of  target  velocities.  Although 
the  initial  estimate  of  target  velocities  is  usually  fairly  coarse, 
it  is  generally  in  the  right  direction  and  results  in  good  estimates 
after  3  to  5  iterations.  Subsequently,  the  target  is  tracked 
very  well  with  only  a  nominal  amount  of  slew.  More  importantly, 
the  computational  speed  is  such  that  it  is  feasible  for  real-time 
implementation,  with  calculations  naving  been  done  at  25  to  100 
iterations  per  second  on  various  computers,  including  time  spent 
in  simulation  support. 

The  most  notable  failure  is  a  high  degree  of  instability  en¬ 
countered  in  dealing  with  real  data  in  the  form  of  digitized 
images.  The  available  image  data,  however,  did  not  have  a  suit¬ 
able  dynamic  range  in  comparison  to  the  noise  level.  Considerable 
improvement  resulted  by  expanding  the  contract  and  filtering  to 
obtain  a  greater  dynamic  range. 

The  results  of  performing  the  above  algorithm  on  simulated 
data  is  presented  in  appendix  A. 


III.  Alternate  Formulation  via  2 


forms . 


The  major  source  of  error  in  the  calculation  of  generalized 
velocities  would  appear  to  be  that  introduced  in  the  numerical 
approximation  of  spatial  derivatives.  Although  the  situation  is 
improved  somewhat  by  filtering  and  the  use  of  multi-point  formulas, 
it  is  still  desirable  to  seek  alternate  approaches.  As  can  be 
seen  from  examination  of  the  algorithm  of  the  preceeding  section, 
any  method  for  calculation  of  generalized  velocities  may  easily 
be  inserted  in  the  basic  tracker. 

In  this  section  we  appeal  to  a  form  of  Stoke 's  Theorem  [1] 
to  obtain  an  integration  based  analogue  of  Theorem  1.3.  The 
formula  obtained  is  strictly  valid  only  when  the  generalized 
velocities  are  constant,  although  is  is  a  useful  approximation 
when  the  rates  of  change  of  the  velocities  are  small. 

We  consider  the  three  dimensional  space  R3  consisting  of 
time  t  and  two  spatial  variables  x  and  y.  Coordinates  £=(t,x,y) 
are  chosen  to  make  a  right-hand  coordinate  system,  and  we  observe 
this  orientation  in  defining  differential  forms.  We  state  the 
form  of  Stoke' s  Theorem  required: 

Stoke1 s  Theorem:  Let  fl  be  a  rectangle  in  (t,x,y)  space  R3  and 

let  w  be  a  differentiable  2-form.  Then 

I  a)  =  /dco  (3.1) 

an  ft 

Here  go  is  of  the  form  co  =  cxgdxdy  +  a,aydt  +  o^dtax,  with 
oig,  a^,  differentiable  functions  of  5eR3,  and 

3a.  3a_ 

=  —  +  a3T  +  dtdxdy> 


i 

J 
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Although  the  results  to  be  presented  may  be  generalized  con¬ 
siderably,  our  derivation  and  experimental  results  will  be  given 
only  for  SA(2) .  Thus,  the  appropriate  vector  fields  and  corres¬ 
ponding  infinitesimal  transformations  may  be  obtained  from  (1.8), 
(1.10)  and  Table  2,  with  x  and  y  substituted  in  the  obvious  manner. 
By  analogy  with  (1.9) ,  for  a  constant  velocity  field 
X  =  (X^,  X^,  X ^ ,  X^)  ,  let  us  define  a  vector  valued  map 
4  2  2 

77:  R  xR  +R  by 


As 


usual,  let 


=  l  X  vi(S)  = 
i=l  1 

2  denote  the  components  of~  . 


r  Ai  +  v  -  v 

x2  +  +  x3y 


(3.2) 


Now,  if  s(t,£)  is  the  trace  corresponding  to  the  generalized 
velocity  field  X  and  F(t,£)  =  f(s(t,0)  is  observed  data,  we  may 
express  Equ.  (1.13)  of  Theorem  1.3  as 


_  3_F  _  7j  3F  ~  3£ 

3t  1  3x  2  3y 

We  intend  to  apply  Stoke ' s  Theorem  to  the 
a)  =  F  dxdy  +^F  dydt  +  dtdx 


(3.3) 

2-form  defined  by 

(3.4) 


The  principal  result  is  stated  as 
Theorem  3.1:  In  the  context  above, 


dw 


+ 


F  dtdxdy 


=  2X^F  dtdxdy 


(3.5) 


To  establish  this,  we  calculate  dw,  using  the  fact  that 
dxdydt  =  dydtdx  =  dtdydx  (whereas,  for  example,  observing  orienta¬ 
tion,  dtdydt  =  -dtdxdy) .  We  have 
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3F  *n  3F 
3t  'l  3x 


M  +  F  air  +  ~2  ly'  +  F  a/)  dtdxd^  3y 


application  of  (3.3)  and  then  (3.2)  this  simplifies  to  dw  = 

/  ^ 

V F  +  F  dtdxdy  =  2X^  F  dtdxdy ,  as  desired. 


By  an  application  of  Stoke’ s  Theorem,  and  a  somewhat  tedious 
calculation,  we  immediately  obtain: 

Theorem  (3.2)  .  Let  n  be  a  rectangle  in  defined  by  opposing 
corners  (t^,  x,  ,  Y 2^  '  In  the  context  described  above,  in  part¬ 
icular  with  X  constant,  we  have 


)  X.k.  =  krt 
iii  1  1  0 


(3.6) 


where  kg,  k.^,...,  k^  are  given  in  Table  3 


Table  3.  Coefficients  resulting  from  Stoke’ s  Theorem. 


k  =  -  J  F  dxdy 

an 

k,  =  /  F  dydt 

1  an 

k-  =  /  F  dtdx 

an 

k^  =  /  xF  dydt  +  /  yF  dtdx  -2  /  F  dtdxdy 

J  an  an  n 

k.  =  /  xF  dtdx  -  /  yF  dydt 

an  an 


It  is  important  that  orientation  be  considered  in  the  evalua¬ 
tion  of  the  coefficients  in  Table  3  (see  [1]).  The  sign  conven¬ 
tion  is  such  that  for  a  principal  2-form  (e.g.,  dxdy)  a  positive 
(negative)  sign  prevails  on  a  face  of  the  rectangle  n  provided 
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that  application  of  a  right-hand  rule  points  outward  from  <in- 

x2 

ward  to  )  the  rectangle  SI.  Writing  for  (similarly  for 

t  and  y)  and  assuming  that  t^  <  t^ ,  x^  <  y^  <  y 2'  ^Y  waY  °f 
example  we  have, 

J  r  dxdy  =  /  /  F  ( t_ , x , y)  dxdy  -  J  J  F  ( t,  ,  x  ,  y )  dxdy  , 

30  y  x  y  x 

and 

/  xF  dydt  =  x-  /  /  F{ t,x- ,y) dydt  -x  /  J  F( t,x. ,y) dydt 
30  tv  *  1  t  y  x 

and 

J  xFdtdx  =  /  /  xF (t,x#y-) dtdx  -  /  J  xF ( t , x , y , ) dtdx . 

30  x  t  x  t 

The  remaining  integrals  may  be  expanded  in  a  similar  fashion. 

Observe  that  differences  are  not  entirely  eliminated  from  the 
final  formulas.  However,  the  formulas  are  so  written  to  indicate 
that  the  differences  are  taken  after  integration,  even  though  in 
certain  cases  the  formula  could  be  collapsed  with  a  difference 
taken  before  evaluation  of  the  iterated  integral. 

The  advantage  of  (3.6)  over  (1.13)  as  a  means  of  calculation 
of  the  generalized  velocities  is  achieved  mainly  by  the  filtering 
effect  of  the  surface  and  volume  integrals.  As  a  matter  of  prac¬ 
tice,  several  rectangles  Q^,...,  0^  are  selected.  Each  rectangle 
0  gives  rise  to  an  equation  of  the  form  (3.6)  , 


4 

1 

i=l 


k(e)  X. 
1  x 


=  k 


(e) 


(3.7) 


The  resulting  system  of  m  equations  in  4  unknowns  may  then  be 
solved  by  a  least-squares  method.  Mote  that  this  approach  may  be 
applied  to  the  feedback  tracking  algorithm  presented  in  Section  2 
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to  calculate  the  velocities  y y  ,  of  the  target  relative  to 
the  window,  replacing  the  corresponding  calculations  based  on 
(2.5) .  This  has  been  implemented  in  a  computer  program  and  tested 
on  real  image  data.  The  results  are  very  encouraging  and  are  pre¬ 
sented  in  part  in  Appendix  S.  This  method  involves  more  com¬ 
putational  overhead,  with  the  best  rate  achieved  to  this  point 
being  about  10  iterations  (  =frames)  per  second.  With  some  stream¬ 
lining  we  believe  that  real-time  rates  of  30  frames  per  second 


can  be  achieved. 
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IV-  Summary  and  Conclusions. 

This  report  presents  a  method  based  on  the  theory  of  Lie 
groups  for  velocity  tracking  in  a  dynamic  image  in  which  the 
motion  of  picture  parts  can  be  ascribed  to  affine  transformations. 

A  feedback  tracking  algorithm  was  developed  and  tested  on  sim¬ 
ulated  data. 

Since  the  random  disturbances  in  real  images  preclude  the 
use  of  simple  methods  for  obtaining  equations  involving  the 
velocities  of  trajectory,  a  method  based  on  integration  of 
differential  forms  was  developed.  This  method  was  incorporated  in 
the  feedback  tracker  and  tested  on  real  image  data.  The  results 
are  very  encouraging,  with  computation  speeds  approaching  ten 
frames  per  second  on  a  VAX  11/780.  We  believe  that  this  method 
is  viable  as  a  component  of  a  real-time  video  tracking  system. 

Among  the  problems  left  outstanding,  a  satisfactory  algorithm 
for  target  acquisition  has  not  yet  been  developed.  In  the  ex¬ 
perimental  work  performed,  the  initial  target  location  was  supplied 
as  an  input  parameter.  To  be  useful,  a  method  for  automating  this 
step  is  essential. 

In  addition,  we  continue  to  experience  problems  with  numerical 
percision.  This  seems  to  be  related  to  the  absence  of  sufficient 
dynamic  range  in  real  image  data,  indicating  that  this  could  be 
improved  by  changes  in  the  capability  of  sensors.  We  feel  that  a 
great  improvement  would  result  from  greater  contrast  in  the  image 
data. 

Finally,  the  class  of  motions  considered  herein  (affine  or 
restricted  affine)  is  not  general  enough  for  many  applications. 
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and  the  methods  need  to  be  extended  to  include  projective  dis¬ 
tortion  as  well. 

The  equations  which  relate  generalized  velocities  to  time- 
varying  images  have  other  applications.  They  have  been  applied 
to  a  problem  in  pattern  matching  with  considerable  success,  as 
fully  described  in  [1] .  The  theoretical  results  and  a  summary 
of  the  experimental  results  of  [11]  have  been  submitted  for 
publication  as  reference  [8] ,  which  is  attached  as  Appendix  C. 
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APPENDIX  A 


Feedback  Tracker  Simulation 

A  typical  imaging  system  might  include  a  sensor  with  a  dia¬ 
meter  (or  cross  section)  of  25  mm  and  an  optical  focal  length  of 
200  mm.  At  a  500  x  500  pixel  density  we  obtain  a  conversion  fac¬ 
tor  of  .05  iran/pix  or  20  pix/mm.  With  a  target  range  of  1  km,  say, 
then  we  obtain  a  conversion  rate  from  sensor  to  target  of  5  m/mm 
3  1  km. 

In  the  results  to  be  presented,  translation  velocities  may 
be  regarded  as  being  given  in  mm/sec.  Conversion  to  pix/sec  or 
m/sec  at  the  target  may  be  done  by  multiplication  by  the  appro¬ 
priate  factor.  Thus,  a  translation  velocity  of  7  mm/sec  at  the 
sensor  corresponds  to  20  pix/sec  or  to  5  m/sec  at  a  target  having 
a  range  of  1  km.  A  magnification  velocity  of  1/sec  translates  by 
the  same  factor  and  would  therefore  represent  a  velocity  of  5m/sec 
3  1  km  toward  the  sensor.  On  the  other  hand,  rotational  velocity 
may  be  considered  as  given  in  radius/sec. 

In  Table  A-l  we  present  the  output  of  the  tracking  simulator 
(the  output  routines  were  modified  slightly  for  ease  of  presenta¬ 
tion)  .  Note  that  the  target  velocities  (10,  -10,  10,  10)  corres¬ 
pond  to  a  3-D  object  with  a  translational  velocity  of  86.6  m/sec 
3  1  km  (about  194  miles/hour)  which  is  rotating  about  3  revolu¬ 
tions  per  second  about  bore  sight.  The  time  base  was  chosen  as 
100  frames/sec.  Inspection  of  the  last  four  columns  of  Table  A-l 
shows  that  the  target  velocities  have  been  acquired  satisfactorily 
after  only  3  frames,  at  t-.03,  and  subsequently  refined  to  exact 


values . 
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The  tracking  simulation  program,  whose  listing  follows,  is 
capable  of  a  real-time  rate  of  33  frames/sec  on  a  VAX  11/780,  in¬ 
cluding  the  time  spent  simulating  target  motion. 


A4 


100 

cmumutttmmmmmiimtumtiimmttmtiitmiim 

200 

c 

c 

300 

c 

ProSram  Title:  TRACK4.F0R 

400 

c 

500 

c 

Function:  Demonstrate  feedback,  tracker  usina  synthetic 

600 

c 

data.  Provide  simulated  motion  consistina  of 

•^oo 

c 

translation*  maanif ication  and  rotation. 

800 

c 

000 

c 

FroSram  Author!  Thomas  6.  Newman 

1000 

c 

Department  of  Mathematics 

1100 

c 

Texas  Tech  University 

1200 

c 

Lubbock*  Texas  79409 

1300 

c 

1400 

c 

Notice:  Permission  is  herewith  aranted  for  use  of  these 

1500 

c 

proarsms*  in  whole  or  in  part*  for  other  than 

1600 

c 

persona’  or  corporate  Sain. 

1700 

c 

1300 

oumitmntnmumummummmmttmtummmm 

1900 

c 

2000 

COMMON  /WINDOW/  I WSIZE *  WIND ( 5 , 5 ) . SWIND <5 *5 ) *  COORD ( 2  *  5 *5 ) 

2100 

COMMON  /EQU/  NEQU.NUNK.A(9,9) »XLAMDA<9> *B(9) 

2200 

COMMON  /FARMS/  TIME , XYSTEP* TSTEP » FEED< 6) 

2300 

COMMON  /C0CHG3/  XTRANI * YTRANI » ECDSI *ESINI » 

2400 

1 

XTRANW* YTRANW*EC0SW*ES1NW* 

2500 

XVELI * YVELI *VHAGI *VR0TI* 

2600 

3 

XVELU*  YVELW*VMAGW* VR0TW 

2700 

100 

TYPE  1. 'ENTER  VELOCITIES  AND  PRESS  RETURN' 

2300 

READ  (5»**END=200>  XVELI * YVELI *VMAGI *VR0TI 

2900 

CALL  I N I T 

3000 

00  175  1=1*15 

3100 

CALL  SAMPLE 

3200 

CALL  DERIV 

3300 

CALL  LINEQ 

3400 

CALL  UPDATE 

3500 

CALL  MOVER 

3600 

CALL  COMPAR 

3700 

173 

CONTINUE 

3800 

GO  TO  100 

3900 

200 

STOP 

4000 

END 

4100 

c 

'200 

c 

»»»»»  SUBROUTINE  INIT  «««<< 

4300 

c 

4400 

c 

This 

subroutine  performs  various  initialisation  functions. 

4500 

c 

4600 

c 

4700 

SUBROUTINE  INIT 

4800 

COMMON  /WINDOW/  IWSIZE»WIND(5»5) *SWIND(5*5) *C00RD(2*5*5) 

4900 

COMMON  /EQU/  NEQU*NUNK*A(9*9) *XLAMDA(9) *B(9) 

5000 

COMMON  /FARMS/  TIME*XYSTEP*TSTEP*FEED<6> 

5100 

COMMON  /COCHGS/  XTRANI *YTRANI *EC0SI *ES1NI * 

3200 

1 

XTRANW* YTRANW*ECOSW *  ESI NW* 

5300 

XVELI *YVELI>VMAGI*VROTI* 

5400 

3 

XVELW  *  YVELW  *  VMAGW*  VROTU 

^  in —  - trliii  n— 


A5 


5500 

C 

5o00 

C 

Da 1 a  is  to  be  sampled  at  each  point  of  a 

5700 

c 

window  of  size  IUSIZE  by  IUSIZE.  A  linear  eouation 

5800 

c 

is  to  be  formed  at  each  ‘interior*  point?  using  the 

5900 

c 

boundary  points  only  in  the  dif ferentiation  process. 

6000 

c 

6100 

c 

The  number  of  eauations?  NEQU?  is  therfore  the  number 

6200 

c 

of  interior  points?  and  involve  NUNK  unknowns?  =  4 

6300 

c 

in  this  program?  but  changeable  to  6  for  GA<2). 

6400 

c 

6500 

c 

XYSTEP  and  TSTEP  are  logical  steps  in  space  and  time. 

6600 

c 

6700 

IUSIZE  =  5 

6800 

NEQU  =  (IUSIZE  -  2)*t2 

6900 

HUNK  =  4 

7000 

XYSTEP  =  0.05 

7100 

TSTEP  =0.01 

7200 

c 

7300 

c 

7400 

c 

GENERATE  THE  COORDINATES  OF  THE  SAMPLE  GRID  RELATIVE  TO 

7500 

c 

THE  WINDOW. 

7600 

c 

7700 

MID  =  IUSIZE/2  +  1 

7800 

DO  10  1=1 ? IUSIZE 

7900 

DO  15  J=1 ? IUSIZE 

3000 

c 

8100 

c 

WINDOW  ORIGIN  AT  CENTER  OF  SAMPLE  GRID 

8200 

c 

3300 

COORD ( 1 ? I ? J)=XYSTEP*FLOAT ( I  -  MID) 

8400 

C00RD(2? I ? J)=XYSTEP*FLOAT ( J  -  MID) 

8500 

c 

8600 

c 

WINDOW  ORIGIN  AT  CORNER  OF  SAMPLE  GRID 

8700 

c 

8800 

c 

COORIK  1 ? I ?  J ) =XYSTEP*FLOAT ( I ) 

8900 

c 

COORD <  2?I?J)=XYSTEP#FL0AT(J) 

9000 

c 

9100 

15 

CONTINUE 

9200 

10 

CONTINUE 

9300 

c 

9400 

c 

SET  FEEDBACKS  TO  UNITY.  LOWER  VALUES  PRODUCE 

9500 

c 

SLOWER  ACQUISITION  BUT  GREATER  STABILITY.  LARGER  VALUES 

9600 

c 

MAY  SPEED  AQUISITION  BUT  CAUSE  INSTABILITY. 

9700 

c 

9800 

FEED(l)  =  1.0 

9900 

FEED(2)  =  1.0 

10000 

FEEIK3)  =  1.0 

1C100 

FEED(4)  =  1.0 

10200 

FEEB(5)  =  1.0 

10300 

FEED ( 6)  =  1.0 

10400 

c 

10500 

c 

INI  TALI ZE  THE  TRAJECTORY  OF  THE  IMAGE  AT  BORE  SIGHT 

10600 

c 

10700 

XTRANI  =  0.0 

10300 

YTRANI  =  0.0 

♦ 


J 
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10900 

ECOSI  =  1.0 

11000 

ESINI  =  0.0 

11100 

C 

11200 

C 

INITALIZE  THE  POSITION  OF  THE  WINDOW  AT  BORE  SIGHT 

11300 

C 

11400 

XTRANW  =0.0 

11500 

YTRANW  =  0.0 

11600 

ECOSW  =1.0 

11700 

ESINW  =  0.0 

11800 

C 

11900 

C 

INITALIZE  THE  WINDOW  VELOCITY  TO  XERO 

12000 

c 

12100 

XVELW  =0.0 

12200 

YVELU  =  0.0 

12300 

VMAGW  =  0.0 

12400 

VROTW  =0.0 

12500 

c 

12600 

c 

GET  AN  INITIAL  SAMPLE  FROM  WINDOW 

12700 

c 

12800 

TIME  =  0.0 

12900 

CALL  SAMPLE 

13000 

c 

13100 

c 

TAKE  THE  INITIAL  STEP  IN  TIME 

13200 

c 

13300 

CALL  MOVER 

13400 

c 

13500 

c 

PRINT  PAGE  HEADINGS 

13600 

c 

13700 

WRITE  ( 6 * 1000 ) 

13300 

WRITE  (6*1010) 

13900 

1000 

FORMAT ( ' 1 ' . 'TIME' *5X. 'XTRANI' *8X. 'YTRANI ' r8X» 'ECOSI '  *9X* 

14000 

1 

' ESINI '  *9X» 'XVELI ' *9X  * ' YVELI '  *9X  * ' VMAGI ' *9X* ' VROTI 

14100 

1010 

FORMAT< 15Xf 'XTRANW' *8X. 'YTRANW' *8X. 'ECOSW' *9X. 

14200 

1 

'ESINW' *9X* 'XVELW' *9X* ' YVELW' *9X* 'VMAGW' *9X* 'VROTW' 

14300 

c 

14400 

c 

PRINT  THE  INITIAL  COMPARISON  BETWEEN  TRAAJECTQRIES 

14500 

c 

14600 

CALL  COMPAR 

14700 

RETURN 

14800 

END 

14900 

c 

15000 

c 

»»»»>  SUBROUTINE  SAMPLE  <«< 

15100 

c 

15200 

c 

This  subroutine  generates  values  in  a  rectangular 

15300 

c 

grid  in  the  tracking  window*  saving  the  old  values. 

15400 

c 

15500 

SUBROUTINE  SAMPLE 

15600 

COMMON  /WINDOW/  IWSIZE»WIND<5*5) *SWIND(5*5) *C00RD(2»5*5 

15700 

DO  20  1=1* IWSIZE 

15800 

DO  25  J=1 * IWSIZE 

15900 

X  =  COQRB< 1  *  I  *  J  > 

loOOO 

Y  =  C00RD(2* I  *  J) 

16100 

SWIND< I  * J)  =  WIND! I  *  J ) 

16200 

HINDU. J)  =  FWIND(X*  Y) 

16300 

16400 

16500 

16600 

16700 

16800 

16900 

17000 

17100 

17200 

17300 

17400 

17500 

17600 

17700 

17800 

17900 

18000 

18100 

13200 

18300 

13400 

13500 

18600 

18700 

13800 

13900 

19000 

19100 

19200 

19300 

19400 

19500 

19600 

19700 

19300 

19900 

20000 

20100 

20200 

20300 

20400 

20500 

20600 

20700 

20800 

20900 

21000 

21100 

21200 

21300 

21400 

21300 

2 1  o  0  0 


CONTINUE 

CONTINUE 

RETURN 

END 


»»»»»  SUBROUTINE  DERIV  ««««< 


This  routine  calculates  the  various  derivatives  needed 
for  formation  of  the  linear  systea  for  the  aeneralized 
velocities  . 

SUBROUTINE  DERIV 

COMMON  /WINDOW/  IUSIZE . WIND ( 5 . 5 ) . SWIND ( 5 . 5 ) . COORD ( 2 . 5 . 5 ) 
COMMON  /EOU/  NEQU»NUNK.A(9.9) »XLAMDA<9) .  B(9> 

COMMON  /FARMS/  TIME . XYSTEP > TSTEP , FEED < 6 ) 

SCALER  =  2.0  *  XYSTEP/TSTEP 
K  =  IUSIZE  -  2 
DO  20  1=1. K 
II  =  I  t  1 
DO  10  J=1,K 
JJ  =  J  M 
L  =  ( I  —  1 )  *K  +  J 
X  =  rGORDd.II.JJ) 

Y  =  C00RD( 2. 1 1 . JJ) 

DX=WIND ( 1 1  +  1  *  J J)  -  WIND(II-l.JJ) 

DY=UIND< I I . J J+l )  -  WIND(II.JJ-l) 

A(L. 1 )=  DX 
A  ( L . 2 ) =  DY 
A ( L . 3 ) =  X*DX  +  Y*DY 
A  (L » 4 )  =  X*DY  -  Y*DX 

B(L)=  SCALER  *  (WIND(II.JJ)  -  SWIND(II . JJ) ) 

CONTINUE 

CONTINUE 

RETURN 

END 

»»»»»  SUBROUTINE  LINEQ  ««<«« 

Modified  froa  the  arauaent  fora! 

LINEQ (M.N.A.X.B.CC) 

SUBROUTINE  LINEQ 

COMMON  /EQU/  NEQU-NUNK. A(9.9) . X < 9 ) . B < 9 > 

INTEGER  CC 

SOLVE  AX=B.  T  HOLDS  AN  UPPER  TRIANGULAR  MATRIX  WHILE  S 
IS  WORKSPACE.  THE  METHOD  FACTORS  A=U«T  WHERE  THE  COLUMNS  OF 
U  ARE  ORTHOGANAL  AND  T  IS  TRIANGULAR.  THE  RESULTING  SYSTEM 
T*X=B '  IS  EASILY  SOLVED  BY  SACK  SUBSTITUTION.  ASSUME  M 
EQUATIONS  AND  N  UNKNOWNS.  (  N  '.  =  N  <  =  9  ) 

THE  MATRIX  OF  COEFFICIENTS.  A  IS  STORED  IN  THE  FIRST  N  ROWS 
AND  THE  FIRST  M  COLUMNS  OF  THE  9X9  A  ARRAY.  THE  ROUTINE 
BRINGS  IN  THE  WHOLE  9X9.  BUT  ONLY  USES  Ad.l)  TO  A(N.M) 

(RECALL  THAT  FORTRAN  STORES  THE  ARRAY  COLUMN-WISE.  BUT 
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27100  C 


27200  C 
27300  C 
27400  C 
27500  C 
27600  C 
27700  C 
27800  C 
27900 
28000 
23100 
28200 
23300 


»»»»»  SUBROUTINE  UPDATE 


«««< 


This  routine  updates  the  velocities  of  the  window 
following  the  calculation  of  the  target  velocities 
relative  to  the  window.  Sensitivita  maa  be  varied 
ba  the  feedback,  factors  in  the  arraa  FEED. 

SUBROUTINE  UPDATE 

COMMON  /UINDOW/  IWSIZEfWIND(5f5> fSWIND(5f5> fC00RD(2f5f5> 
COMMON  /EOU/  NEQU  fNUNK  f A<  9  >  9  ) t XLAMDA( 9 ) r  B <  9  > 

COMMON  /FARMS/  TIME f XYSTEP f TSTEP f FEED ( 6 ) 

COMMON  /COCHGS/  XTRANI f YTRANI fECOSI fESINI f 


28400 

1 

XTRANWfYTRANUfECOSUfESINWf 

28500 

n 

XVELIfYVELIfVMAGIfVROTIf 

23600 

3 

XVELW f  YVELW fVMAGU f  VROTU 

23700 

DVELX  =  (ECOSW  *  XLAMDA(l)  -  ESINW  *  XLAMDA( 2) )  - 

29800 

1 

(XLAMDA<3)  *  XTRANW  -  XLAMDAC4)  *  YTRANU) 

28900 

DVELY  =  (ESINW  *  XLAMDA(l)  P  ECOSW  *  XLAMDA(2)>  - 

29000 

0 

(XLAMDA(4>  t  XTRANU  +  XLAMBA(3)  *  YTRANW) 

29100 

DMAGV  =  XLAMDA<3) 

29200 

DVROT  =  XLAMDAC4) 

29300 

C 

29400 

XVELW  =  XVELW  -  FEED(l)  *  DVELX 

29500 

YVELW  =  YVELW  -  FEED<2)  *  DVELY 

29600 

VMAGW  =  VMAGW  -  FEEDC3)  *  DMAGV 

29700 

VRQTW  =  VROTW  -  FEED(4)  *  DVROT 

29800 

C 

29900 

RETURN 

30000 

END 

30100 

C 

30200 

C 

»»»»»  SUBROUTINE  MOVER  «««« 

30300 

c 

30400 

c 

This  routine  moves  the  window  ba  taking  a  step  in 

30500 

c 

the  differential  eauations  for  the  affine  tranformation 

30600 

c 

which  controls  the  window  location.  The  method  used  is 

30700 

c 

a  simple  Euler  method. 

30800 

c 

30900 

c 

The  routine  also  simulates  the  motion  of  the  target  ba 

31000 

c 

solving  the  corresponding  differential  eauation  for  the 

31100 

c 

target.  This  portion  would  be  removed  if  real  data  were 

31200 

c 

being  used. 

31300 

c 

31400 

31500 

31600 

31700 

21300 

31900 


SUBROUTINE  MOVER 

COMMON  /FARMS/  TIME f XYSTEP f TSTEP t FEED ( 6 ) 
COMMON  /COCHGS/  X TRAN  I f YTRANI fECOSI fESINI » 

1  XTRANWfYTRANWfECOSWfESINWf 

2  XVELIfYVELIfVMAGIfVROTIf 

3  XVELW f  YVELW f VMAGW f VROTW 


32000 

DECOS 

=  VMAGW 

* 

ECOSW 

-  VROTW 

*  ESINW 

32100 

DE3IN 

=  VROTU 

* 

ECOSW 

+  VMAGW 

*  ESINW 

32200 

ECOSW 

=  ECOSW 

+ 

TSTEP 

*  I'ECOS 

32300 

32400  C 

E3INU 

=  ESINW 

TSTEP 

*  DE3IN 

32500 

DXTRAN  =  XVELU  +  VMAGU4XTRANW  -  VROTU*YTRANU 

32600 

DYTRAN  =  YVELU  +  VRQTWtXTRANW  4  VMAGU*YTRANW 

32700 

XTRANU  =  XTRANW  +  DXTRANITSTEP 

32300 

YTRANW  =  YTRANW  t  DYTRAN*TSTEP 

32900 

C 

33000 

C 

The  following  portions  simulate  motion  of  the  target. 

33100 

c 

33200 

DECOS  =  VMAGI  *  ECOSI  -  VROTI  *  ESINI 

33300 

DESIN  =  VROTI  *  ECOSI  +  VMAGI  *  ESINI 

33400 

ECOSI  =  ECOSI  +  TSTEP  *  DECOS 

33500 

ESINI  =  ESINI  +  TSTEP  *  DESIN 

33600 

c 

33700 

DXTRAN  =  XVELI  +  VMAGItXTRANI  -  VROTUYTRANI 

33800 

DYTRAN  =  YVELI  +  VROTI*XTRANI  +  VMAGIIYTRANI 

33900 

XTRANI  =  XTRANI  +  DXTRANCTSTEP 

34000 

YTRANI  =  YTRANI  +  DYTRANCTSTEP 

34100 

c 

34200 

c 

Increment  time. 

34300 

c 

34400 

TIME  =  TIME  +  TSTEP 

34500 

c 

34600 

RETURN 

34700 

END 

34800 

c 

34900 

c 

»»»»»  SUBROUTINE  COMPAR  ««««« 

35000 

c 

35100 

c 

This  routine  produces  printed  output  for  evaluation 

35200 

c 

purposes?  and  is  therefore  ancillary  to  the  operation 

35300 

c 

of  the  tracker. 

35400 

c 

35500 

SUBROUTINE  COMPAR 

35600 

COMMON  /PARMS/  TIMErXYSTEP? TSTEP .FEED ( 6 > 

35700 

COMMON  /COCHGS/  XTRANI ? YTRANI »EC0SI ? ESINI » 

35800 

1 

XTRANU? YTRANW rECOSW.ESINW? 

35900 

2 

XVELI ? YVELI .VMAGI ? VROTI * 

36000 

3 

XVELU?  YVELU .VMAGWrVROTW 

36100 

WRITE (6. 2000) TIME. XTRANI .YTRANI .ECOSI .ESINI .XVELI .YVELI .VMAGI? 

36200 

1 

VROTI 

36300 

2000 

FORMAT (2X?F3.2?8(3X?E11,4)) 

36400 

WRITE! 6? 2010) XTRANU? YTRANW? ECOSW?ESINW?XVELW? YVELW.VMAGU? 

36500 

1 

VROTW 

36600 

2010 

F0RMAT(10X?8(3X?E11.4)) 

36700 

RETURN 

36800 

END 

36900 

C 

37000 

C 

>>'»»>»  FUNCTION  FWIND  <<<<<<<<<< 

37100 

c 

37200 

c 

This  function  returns  a  value  at  the  point  (x.y) 

37300 

c 

in  the  window.  It  first  maps  to  the  corresponding 

37400 

c 

point  in  absolute  image  coordinates  and  calls  for 

37500 

c 

image  value  at  that  point  (see  FIMAGE  below). 

37o00 

37700 

37300 

c 

FUNCTION  FVI NB( X  ?  Y  ) 

COMMON  /COCHGS/  XTRANI -YTRANI .ECOSI  .ESINI ? 

1 


FUNCTION  FVINB<X?Y> 

COMMON  /COCHGS/  XTRANI . YTRANI fECOSI .ESINI . 


All 


37900 

1 

XTRANW*  YTRANW*ECOSW*ESINW* 

39000 

n 

XVELI* YVELI *UMAGI *URQTI * 

3S100 

3 

XVELU  *  YVELW  * VMAGW * URQTW 

33200 

c 

38300 

c 

38400 

c 

33500 

XlhAGE  =  XTRANW  +  £COSW*X  -  ESINW*Y 

38600 

YIMAGE  =  YTRANW  +  ESINW*X  +  ECOSU*Y 

33700 

FWIND  =  F IMAGE (X  IMAGE f YIMAGE) 

38800 

RETURN 

38900 

END 

39000 

c 

39100 

c 

»»»»»  FUNCTION  FIMAGE  ««««« 

39200 

c 

39300 

c 

This  function  returns  the  value  at  point  (xiv)  in  the 

39400 

c 

absolute  image  coordinate  system.  It  maps  to  target 

39500 

c 

coordinates  and  calls  for  the  value  at  the  correseondins 

39600 

c 

point  on  the  tarset  (see  FQBJ  below). 

39700 

c 

39800 

FUNCTION  FIMAGE(X* Y) 

39900 

COMMON  /COCHGS/  XTRANI » YTRANI »EC0SI *ESINI * 

40000 

1 

XTRANW* YTRANU*£COSU* ESI NU* 

40100 

XUELI » YUELI *VMAGI » VROTI * 

40200 

3 

XVELW* YVELW*UMAGU* VROTW 

40300 

c 

40400 

c 

40500 

DET  =  EC0SI**2  +  ESINI**2 

40600 

DX  =  X  -  XTRANI 

40700 

BY  =  Y  -  YTRANI 

40800 

c 

40900 

c 

41000 

XOB  J  =  (  EC0SI*DX  +  ESINItDY) /DET 

41100 

YOB J=  (  -ESINUDX  4  EC0SI*DY) /DET 

41200 

FIMAGE  =  FOBJ(XOBJ.YOBJ) 

41300 

RETURN 

41400 

END 

41500 

c 

41o00 

c 

»»»»»  FUNCTION  FOBJ  <<<<<<<<<< 

41700 

c 

41800 

c 

This  function  returns  the  gray  value  at  a  point  (x*y> 

41900 

c 

in  target  coordinates.  For  actual  tracking*  this 

42000 

c 

routine  would  be  replaced  by  one  which  retrieves  from 

42100 

c 

the  image  database.  In  the  simulation*  however*  the 

42200 

c 

routine  merely  returns  a  synthetic  value  generated 

42300 

c 

from  an  expression. 

42400 

c 

42500 

FUNCTION  FOBJ(X.Y) 

42600 

FOBJ  =  1.0  t  10.0*X  -5 . 0*Y  t  20.0*X*Y 

42700 

RETURN 

42800 

END 

* 


APPENDIX  3 


Tracking  With  Differential  Forms 

A  tracking  program  was  developed  which  utilized  the  theory 
presented  in  Section  III.  In  order  to  test  this  program,  digit¬ 
ized  video  images  were  obtained  from  the  Advanced  Technology 
Office,  Instrumentation  Directorate,  White  Sands  Missile  Range. 

One  such  image  is  shown  in  Figure  B-l.  A  sequence  of  test  images 
was  prepared  from  this  data  by  shifting  to  inject  additional  mo¬ 
tion.  Sections  of  the  first  six  frames  are  shown  in  the  left  hand 
column  of  Figure  B-2. 

Parameters  in  the  tracking  program  were  set  to  use  a  3x3  win¬ 
dow  in  3  consecutive  frames  to  form  a  3x3x3  rectangle  in  space. 
Since  the  pr'  ;ram  uses  9  such  windows,  the  actual  track  gate  con¬ 
sisted  of  9x9x3  points. 

The  right  hand  column  of  Figure  B-2  shows  the  output  frames 
obtained  by  the  tracker.  Observe  that  the  output  lags  the  input 
by  the  depth  of  the  track  gate.  This  is  why  there  are  fewer  out¬ 
puts  frames  than  input  frames.  We  see  that  the  target  was  ac¬ 
quired  immediately,  and  successfully  tracked  over  the  full  se¬ 
quence  of  frames. 

The  computation  rate  was  about  10  frames  per  second  on  a 
VAX  11/780.  However,  the  program  is  written  to  allow  selection 
of  window  sizes  and  could  be  streamlined  a  great  deal. 

Although  the  results  shown  in  Figure  3-2  are  impressive,  we 
hasten  to  point  out  that  the  motion  is  mostly  translation,  and  our 
attempts  to  test  the  algorithm  on  a  wide  range  of  motions  have 
been  frustrated  by  availability  of  data.  Further  work  is  ongoing, 
and  refinements  to  the  tracking  program  are  expected  to  be  forth¬ 


coming. 


B4 


20"  C  PROGRAMMER  !  DONNA  K.  TERRAL 

40  C  DEPARTMENT  OF  MATHEMATICS 

60  C  TEXAS  TECH  UNIVERSITY 

30  C 

100  C  PERMISSION  IS  HEREWITH  GRANTED  TO  UTILIZE  THIS  PROGRAM  FOR 

120  C  OTHER  THAN  PERSONAL  OR  CORPORATIVE  GAIN. 

140  C 


160  C 

130  C  ******************************************************************** 

200  C  ***************************  MAIN  PROGRAM  *************************** 

220  C  ******************************************************************** 

240  C  PURPOSE!  GIVEN  A  SET  OF  IMAGES*  TRACK  A  TRAGET  BY  USING 

260  C  INTEGRATION  TO  DETERMINE  THE  MOTIONS  ( TRANSLATION* ROT AT I  ON* 

280  C  MAGNIFICATION)  OF  THE  TARGET. 

300  C  ******************************************************************** 

320  C 


340  C 
360  C 
330 

400  C 
420 
440 
460 
430 
500 
520 
540 
560 
530 
600 
620 
640 
660 
630 
700 
720 
740 
760 

780  C 
300  C 
320 
340 
360 

330  C 
900 

920  C 
’40  C 
360 

930  C 
1000  C 


*  COMMON  DECLARATIONS  * 

PARAMETER  ( IU=9 . JW=9 . KW=3» IB=3. JB=3»KB=3. ISIZE=128 . JSIZE=128) 

COMMON  /WINDOW/  WIND ( IW* JW .KW) . CORRD ( ISIZE . JSIZE » 2) 

INTEGER  WIND 

COMMON  /ORIGIN/  10 , JO . KO . XO .  YO 
INTEGER  10 » JO » KO .  XO  >  YO 
COMMON  /BUFFER/  BUF ( IB . JB.KB > 

INTEGER  BUF 

COMMON  /CORNER/  XI , X2 * Y1 , Y2 »T1 . T2 
INTEGER  XI >  X2 .  Y1 . Y2 . T1 » T2 

COMMON  /WEIGHT/  W(  IB .  JB » KB  > » WX  ( JB»  KB) » WY !  IB  .KB) » WT  ( IB .  JB ) 
INTEGER  W.WX.WY.WT 

COMMON  /COCHGS/  XVELI * YVELI .VROTI .VMAGI .TIME.TSTEP 
REAL  XVELI . YVELI . VROTI .VMAGI . TIME . TSTEP 
COMMON  /PARMS/  ROW. COL. PROW. PCOL.NUMROW.NUMCOL.WTEST 
INTEGER  ROW. COL. PROW. PCOL.NUMROW.NUMCOL.WTEST 
COMMON  /IMAGE/  IMAGE! ISIZE. JSIZE) 

INTEGER*2  IMAGE 

COMMON  /EQU/  ALPHA ( 4 > . BETA.COEFF ( 9 . 4 ) . VECTOR ( 9 > , SOLN ( 4 ) . FEED ( 4 ) 
INTEGER  ALPHA. BETA 

*  MAIN  VARIABLES  * 

INTEGER  I1.J1.K1 
INTEGER  CC. COUNT 
DATA  COUNT. CC.K1  /l.l.l/ 

CALL  INIT 

***  MAIN  LOOP  **** 

DO  100  LOOP  =1.4 

ALL  OR  A  PORTION  OF  THE  INTEGRATION  RESULTS  CAN  BE  USED 


1020  SOLN ( 1 )  =  FEED!  1 )  *  SOLN(l) 

1040  SOLN! 2)  =  FEED! 2)  *  S0LN(2) 


1060  SOLN! 3)  =  FEED!3)  *  S0LN(3> 

1080  SOLN! 4 )  =  FEED! 4 )  *  S0LN<4> 
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CALL  WINDOW  (LOOP) 

CALL  PRINT 

***  LOOP  COMPUTES  MOTION  USED  TO  MOVE  WINDOW  *** 
DO  1  J1  =  l.JUfJB 
DO  2  II  =  1 r IW» IB 

CALL  GETBUF  (I1.J1,K1> 

IF  (LOOP  .NE.  i>  GO  TO  10 
CALL  BLOW 
CALL  GETEQU 
CALL  MATRIX  (COUNT) 

COUNT  =  COUNT  +  1 
CONTINUE 
CONTINUE 

CALL  LINED  (COEFF»SOLN,VECTOR»9,4»CC> 

COUNT  =  1 

tttt  END  MOTION  LOOP  ttt 

TIME  =  TIME  +  TSTEP 
CONTINUE 

***  END  MAIN  LOOP  t*« 


C  tt****tt*t**t*t*tttt*tt*t*ttt%tt*t*t*tf.ttt**ttt***t*ttt*tttttttttt* 

C  SUBROUTINE  GETBUF  FILLS  A  BUFFER  ARRAY  WHICH  WILL  CONTI vN  THE 

C  DATA  POINTS  FORMING  THE  CUBE  TO  BE  USED  IN  SUBROUTINE  GETEQU. 

C  THE  POSITION  IN  WHICH  TO  BEGIN  IS  PASSED  THRU  THE  ARGUMENTS. 

C  <BUF ( 1 » 1 » 1 ) =WIND ( 1 1 r  J1 f  K 1 ) > 

c  mmmmmxx***mxxmxxm**m*xxx***xxxxtm*t*xxt**x**tt* 

c 

SUBROUTINE  GETBUF  <I1,J1»K1) 

C 

C  *  ARGUMENTS  * 

INTEGER  1 1 r J1 r K1 
r 

C  *  COMMON  DECLARATIONS  * 

C 

PARAMETER  ( IW=9 » JW=? , KW=3 r I B  =  3f JB=3 r KB=3» ISIZE=128 * JSIZE= 128 ) 


COMMON  /WINDOW/  WIND ( IW, JW.KW) ,CORRD( ISIZE » JSIZE. 2 ) 
INTEGER  WIND 

COMMON  /ORIGIN/  10 . JO » KO f XO » YO 
INTEGER  IOr JO»KO»XO»YO 
COMMON  /BUFFER/  BUF( IB, JBrKB) 

INTEGER  BUF 

COMMON  /CORNER/  XI »X2» Y1 .Y2fT1 »T2 
INTEGER  XI »X2» Y1 r Y2 »T1 » T2 


*  LOCAL  VARAIBLES  * 


r 

i 

B6 

2130 

INTEGER  X  *  Y  *  T 

2200 

C 

2220 

XI  =  11  -  10 

2240 

o 

1 

r4 

II 

•H 

>- 

2260 

T 1  =  K 1  —  K0 

2230 

X2  =  XI  +  IB  -  1 

;  2300 

Y2=Y1+JB-1 

!  2320 

T2  =  T1  +  KB  -  1 

|  2340 

c 

!  2360 

BO  30  I  =  1  *  IB 

2380 

DO  30  J  =  ltJB 

2400 

DO  30  K  =  l.KB 

?  2420 

r4 

1 

M 

W 

II 

X 

2440 

Y  »  J1  +  J  -  1 

j  2460 

T  =  K1  +  K  -  1 

j  2430 

BUF ( I  * J*K)  =  W1ND(X»Y*T) 

j  2500 

30  CONTINUE 

2520 

RETURN 

■  2540 

END 

I  2560 

C 

i  2580 

C 

tt%t%*%%tt%***t*********%*t***%***tt**t*t********tt*t**tt*****ttt* 

|  2600 

C 

!  2620 

C 

SUBROUTINE  MATRIX  TAKES  THE  ALPHAS  AND  BETA  FROM  THE  SUBROUTINE 

2640 

C 

GETEQU  AND  PUTS  THEM  IN  THE  FORM  AX-^B  WHERE  LOOPS  OF  GETEGU 

j  2660 

C 

FORM  THE  2-DIMENSIONAL  ARRAY  A  AND  THE  VECTOR  B. 

!  2680 

C 

ONCE  AX=B  IS  FORMED  LINED  IS  USED  TO  SOLVE  FOR  X. 

!  2700 

C 

2720 

C 

IN  MATRIX*  COEFF (9*4)  IS  ARRAY  A  AND  VECTOR!?)  IS  ARRAY  B 

|  2740 

c 

tt%%t%%%*%%%%***tt%%*******%**t*%***%*t*ttttttt**tt****ttt*****t*t 

j  2760 

c 

!  2780 

SUBROUTINE  MATRIX  (COUNT) 

|  2300 

c 

2820 

c 

*  ARGUMENTS  * 

!  2340 

INTEGER  COUNT 

j  2860 

c 

2380 

c 

*  COMMON  DECLARATIONS  * 

2900 

COMMON  /EQU/  ALPHA(4) *  BETA* COEFF (9*4) r VECT0R<9) rS0LN( 4 ) r FEED (4 ) 

2920 

INTEGER  ALPHA*BETA 

2940 

c 

2  5  60 

COEFF (COUNT  *  1 )  =  FLOAT (ALPHA! 1 ) ) 

2980 

COEFF (COUNT  *  2)  =  FLOAT! ALPHA(2) ) 

3000 

COEFF (COUNT  *3)  =  FLOAT! ALPHAC3) ) 

3020 

COEFF ( COUNT »  4 )  =  FLOAT < ALPHA< 4 ) ) 

3040 

VECTOR (COUNT )  =  FLOAT (BETA) 

3060 

RETURN 

3080 

END 

3100 

c 

3120 

c 

%ttt%t%ttt%***tt*t**tt*tt**t****tt***t**ttt******ttt*t**ttt*ttt*tttt 

1  3140 

c 

%rn.%*ttttt*t%t**t*****%*tt**t*%***tt*tt*t****tt**t*t****%tt**ttttt*t 

)  3160 

c 

SUBROUTINE  GETEQU 

,  3130 

c 

COMPUTES  CONSTANTS  ALPHA! 1 )  THRU  ALPHA ( 4 )  AND  BETA  IN  THE 

3200 

c 

EQUATION  ALPHA(l)  *  LAMBDA! 1 )  +  ALPHA( 2 )  *  LAMBDA< 2 )  + 

3220 

c 

AIPHA(3>  *  LAMBDA( 3)  4-  ALPHA(4)  *  LAMBDA! 4 )  =  BETA*  WHERE 

3240 

. -  - 

c 

ALPHA ( 1 )  THRU  ALPHA( 4)  AND  BETA  ARE  FORMED  FROM  SURFACE  AND 
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3260  C  VOLUME  INTEGALS  OVER  A  CUBE  INDEXED  BY  X1fX2fY1fY2fT1fT2. 

3280  C  THESE  INTEGALS  ARE  NUMERICALLY  INTEGRATED  USING  THE  TRAPEZIQD 

3300  C  RULE. 

3320  C 

3340  C  let: 

3360  C  FXYT  REPRESENT  THE  VOLUME  INTEGAL  OVER  THE  CUBE 


3380 

C 

FXY  1 

REPRESENT 

THE 

SURFACE 

INTEGAL 

OVER  THE 

FACE 

XY  9 

T=1 

3400 

C 

FXY2 

■ 

• 

• 

a 

• 

FACE 

XY  9 

T=2 

3420 

C 

FYT1 

1 

• 

■ 

a 

1 

FACE 

YT  9 

X*1 

3440 

c 

FYT2 

« 

• 

> 

• 

a 

FACE 

YT  8 

X=2 

3460 

c 

FTX1 

• 

■ 

■ 

a 

a 

FACE 

TX  9 

Y=1 

3480 

c 

FTX2 

• 

■ 

a 

a 

a 

FACE 

TX  9 

Y=2 

3500 

c 

YFYT1 

REPRESENT 

THE 

INTEGAL 

OVER 

Y 

* 

(FACE 

YT)  9 

X=1 

3520 

c 

YFYT2 

• 

• 

• 

• 

Y 

t 

(FACE 

YT)  9 

X=2 

3540 

c 

XFTX1 

1 

f 

a 

a 

X 

t 

(FACE 

TX)  9 

Y=1 

3560 

c 

XFTX2 

■ 

■ 

a 

a 

X 

* 

(FACE 

TX)  e 

Y=2 

3580  C 

3600  C  THEN*. 

3620  C  ALPHA< 1 )  *  FYT2  -  FYT1 

3640  C  ALPHA<2)  =  FTX2  -  FTX1 

3660  C  ALPHA(3>  =  X2*FYT2  -  XUFYT1  -  2»FXYT  +  Y21FTX2  -  YUFTX1 

3680  C  ALPHAC4)  =  YFYT1  -  YFYT2  +  XFTX2  -  XFTX1 

3700  C  BETA  =  FXY1  -  FXY2 

3720  c  *mm**mtm*mm***t**mm*mm*mmm*mm***m*m 

3740  c 

3760  SUBROUTINE  GETEQU 

3780  C 

3800  C  *  COMMON  DECLARATIONS  * 

3820  C 

3840  PARAMETER  ( IW=9 f JU=9 fKU=3f IB=3 f JB=3fKB=3 f ISIZE=128f JSIZE=128> 

3860  C 

3880  COMMON  /BUFFER/  BUF ( IB » JB »KB ) 

3900  INTEGER  BUF 

3920  COMMON  /CORNER/  XI fX2f Y1 fY2fT1 f T2 

3940  INTEGER  X1fX2fY1fY2fT1fT2 

3960  COMMON  /WEIGHT/  W( IBf JBfKB) f WX( JBfKB) fUY( IBfKB) fWT( IBf JB) 

3980  INTEGER  WfUXfWYfWT 

4000  COMMON  /EQU/  ALPHA<4 ) f  BETAfCOEFF( 9f  4 ) f VECT0R(9) fS0LN(4 )  fFEED( 4) 

4020  INTEGER  ALPHAfBETA 

4040  C 

4060  C  *  LOCAL  VARIABLES  * 

4080  INTEGER  FXY1 fFXY2fFYT1 fFYT2fFTX1 fFTX2f YFYT1 f YFYT2 fXFTXI fXFTX2fFXYT 

4100  C 

4120  FXYT  =  0 

4140  FXY 1  =  0 

4160  FXY2  =  0 

4180  FYT1  =  0 

4200  FYT2  =  0 

4220  FTX1  =  0 

4240  FTX2  =  0 

4260  YFYT1  =  0 

4280  YFYT2  =  0 

4300  XFTX1  *  0 

4320  XFTX2  =  0 
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4340  C 


4360 

4380 

4400 

4420 

4440  50 

4460 

4480 

4500 

4520 

4540 

4560 

4580  60 

4600 

4620 

4640 

4660 

4680 

4700 

4720  70 

4740 

4760 

4780 

4300 

4820  80 

4840  C 

4860 

4880 

4900 

4920 

4940 

4960  C 

4980 

5000 


DO  50  I  =  1» IB 
DO  50  J  =  lfJB 

FXY1  =  FXY1  +  BUF <  I  >  J » 1 )  *  WT(I»J) 
FXY2  =  FXY2  +  BUF(I,J,KB>  *  WT ( I » J ) 
CONTINUE 
DO  60  J  =  1  r  J B 
DO  60  K  =  1  iKB 

FYT1  =  FYT1  +  BUFU,J»K>  *  WX(J,K) 
FYT2  =  FYT2  ¥  BUF(IB»J,K>  *  WX(J,K) 


YFYT1  =  YFYT1  +  CJ+Y1  —  1>  *  BUFC1,J,K)  *  MX(J»K> 
YFYT2  =  YFYT2  +  (J  +  Y1  -  1)  *  BUFUB.JtK)  *  WX(J,K) 
CONTINUE 


DO  70  I  =  1  *  IB 
DO  70  K  =  1»KB 

FTX1  =  FTX1  +  BUF < I » 1 » K )  *  WY<I,K) 

FTX2  =  FTX2  +  BUF(I,JB,K>  *  UY<I»K) 

XFTX1  =  XFTX1  +  (I  +  XI  -  1)  *  BUF ( I > 1  * K)  *  MY(I»K) 
XFTX2  =  XFTX2  +  (I  +  XI  -  1)  *  BUF(I»JB»K)  *  UY(IfK) 
CONTINUE 
DO  80  I  =  If  IB 
DO  80  J  =  l.JB 
DO  80  K  =  1 i KB 

FXYT  =  FXYT  +  BUF(I»J,K>  t  W(I»J,K> 

CONTINUE 


ALPHA! 1 )  =  FYT2  -  FYT1 
ALPHA(2)  =  FTX2  -  FTX1 

ALPHAC3)  =  X2  *  FYT2  -  XI  *  FYT1  -  FXYT  +  Y2  *  FTX2  -  Y1  *  FTX1 
ALPHA(4)  =  YFYT1  -  YFYT2  +  XFTX2  -XFTX1 
BETA  =  FXY 1  -  FXY2 


RETURN 

END 


5040  C 

5060  C  tt**t**t***t*t*tt*********t*ttt$******ttt**tttttt*ttt**ttttttttttttt 


3080  C  SUBROUTINE  BLDW  BUILDS  THE  WEIGHING  ARRAYS  FOR  ’HE 

5100  C  TRAPEZIOD  RULE  USED  IN  THE  SUBROUTINE  GETEQU .  THESE  ARRAYS 

5120  C  ARE  PASSED  TO  GETEQU  THRU  A  COMMON  BLOCK. 

5140  C  t***%*ti(ttt**tt**%%tt%**t%tt*t**%%***t*tt*t*t*t*tt*%*%**t**t*t**%*% 

5160  C 

5180  SUBROUTINE  BLDW 

5200  C 

5220  C  *  COMMON  DECLARATIONS  * 

5240  C 

5260  PARAMETER  ( IW=9, JW=9.KW=3» IB-3* JB=3»KB=3» ISIZE=128» JSIZE=128) 

5230  C 

5300  COMMON  /WEIGHT/  W< IB » JB> KB >  r WX ( JB r KB >  > WY ( IB>KB ) > WT ( IB » JB > 

5320  INTEGER  WfWXrWYrWT 

5340  C 

5360  C  LOCAL  VARAIBLES 

5330  INTEGER  X.YfT 

3400  C 


5420 

X  =  IB  -  1 

5440 

y  =  jb  -  i 

5460 

5430  C 

T  =  KB  -  1 

5500 

DO  40  I  ■  2»X 

5520 

DO  40  J  =  2»Y 

5540 

DO  40  K  3  2.T 

5560 

M(  I > 1 > 1 )  =  2 

5580 

U(I'JB'l)  =  2 

5600 

WUfJrl)  3  2 

5620 

U( IB» J» 1 )  3  2 

5640 

U( 1 • JBrK)  =  2 

5660 

U( IBi 1 >K)  3  2 

5680 

H(IflfKB)  3  2 

5700 

W( 1 » Jf KB>  3  2 

5720 

U(IBrJtKB)  3  2 

5740 

VK I f JBiKB)  3  2 

5760 

U( IB> JBi K)  3  2 

5780 

U( 1 » 1 >K)  3  2 

5800 

U(I» JrK)  3  8 

5820 

U(l.J.K)  3  4 

5840 

W(I.l.K)  3  4 

5860 

W(IfJil)  3  4 

5880 

W< ItJS/K)  -  4 

5900 

U(IB.J.K)  3  4 

5920 

U( I » J>KB)  3  4 

5940 

MY ( I » l )  3  2 

5960 

MY(l.K)  3  2 

5980 

UY(I.K)  3  4 

6000 

MY ( I f KB)  3  2 

6020 

UY(IBtK)  3  2 

"MO 

UT(I»1)  3  2 

,0 

MT(liJ)  3  2 

6080 

MT(IfJ)  3  4 

6100 

MT(IfJB)  3  2 

6120 

UT(IB.J)  3  2 

6140 

MX ( JB  >  K )  3  2 

6160 

MX< JfKB)  3  2 

6180 

MX(lrK)  3  2 

6200 

UX(Jfl)  3  2 

6220 

WXCJ»K)  3  4 

6240  40 

CONTINUE 

6260 

U(l.lrl)  3  1 

6280 

M( 1 > JB  *  1 )  3  1 

6300 

U(liJBfKB)  3  1 

6320 

U ( 1 » 1 >KB  )  3  1 

6340 

W< IB» 1 » 1 )  3  1 

6360 

U(IBtJBfl)  3  1 

6380 

U(IBfI.KB)  3  1 

6400 

W(IB»JBfKB)  3  1 

6420 

MX<1»1>  3  1 

6440 

MX ( 1 »KB)  3  1 

6460 

UX(JBil)  3  1 

6480 

UX(JB»KB)  3  l 

M  w  rj  m  r  J 
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6500  WY ( 1 » 1 )  =  1 

6520  WY( 1 >KB)  =  1 

6540  UY(IB»1)  =  1 

6560  WY< IBfKB)  =  1 

6580  WT ( 1 » 1 )  =  1 

6600  WT<1»JB>  =  1 

6620  WT  < IB » 1 1  =  1 

6640  UT<IB*JB)  =  1 

6660  RETURN 

6680  END 

6700  C 
6720  C 
6740  C 

6760  C  SUBROUTINE  PRINT  WRITES  TO  UNIT  =  66  THE  MOTION  PARAMETERS 

6780  C  OF  THE  IMAGE  AND  THE  WINDOW  AT  EACH  TIME  STEP.  THE  PIXEL 

6800  C  VALUES  IN  THE  WINDOW  MAY  ALSO  BE  PRINTED  IF  NEEPED. 

6820  c 

6840  C 

6860  SUBROUTINE  PRINT 

6880  C 

6900  C  *  COMMON  DECLARATIONS  * 

6920  C 

6940  PARAMETER  ( IW=9» JW=9  ,KW=3.  IB=-3» JB=3»KB=3»ISIZE=128r JSIZE=128> 

6960  C 

6980  COMMON  /WINDOW/  WIND ( IW» JW»KW) * CORRD( ISIZE r JSIZE r 2 > 

7000  INTEGER  WIND 

7020  COMMON  /ORIGIN/  10 » JOi K0>X0>Y0 

7040  INTEGER  10 > JO » KO » XO > YO 

7060  COMMON  /CQCHGS/  XVELI » YVELI » VROTI r VHAGI r TIME jTSTEP 

7030  REAL  XVELI » YVELI » VROTI » VMAGI f TIMEjTSTEP 

7100  COMMON  /PARMS/  ROW,COL»PROWfPCOL>NUMROWpNUMCOL,WTEST 

7120  INTEGER  ROU>COL»PROWfPCOLfNUMROW»NUMCOL»WTEST 

7140  COMMON  /EQU/  ALPHA( 4 ) f BETA, COEFF ( 9 , 4 ) , VECTOR ( 9 ) , SOLN( 4 ) , FEED( 4 ) 

7160  INTEGER  ALPHA, BETA 

7180  C 

7200  C  *  LOCAL  VARIABLES  * 

7220  REAL  XTRANW, YTRANW, ROTW* MAGW, XVELW»YVELW»VRQTW,VMAGW»XTRANI , 

7240  1  MAGI *ROTI , YTRANI 

7260  C 

7280  DATA  XTRANW. YTRANW, ROTW, MAGW  /O, 0,0,0/ 

7300  C 

7320  XTRANI  =  XVELItTIME 

7340  YTRANI  =  YVELI*TIME 

7360  MAGI  =  VMAGItTIME 

7380  ROTI  =  VROTI*TIME 

7400  C 

7420  XTRANW  =  SOLN(l)  +  XTRANW 

7440  YTRANW  =  S0LN<2>  +  YTRANW 

7460  MAGW  =  30LN<3>  +  MAGW 

7480  ROTW  =  SOLN<  4)  +  ROTW 

7500  C 

7520  XVELW  =  SOLN(  1XU/TSTEP) 

7540  YVELW  =  S0LN<2>!< l/TSTEP) 

7560  VMAGW  =  SOLN ( 3 ) I ( l/TSTEP ) 
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7530 

VROTU  =  S0LN(4)*( 1/TSTEP) 

7600 

C 

7620 

IF  (UTEST  .EQ.  1)  THEN 

7640 

ITEMP  =  COL  -  PCOL  +  (IU+D/2 

7660 

JTEMP  *  ROU  -  FROM  +  (JU+l)/2 

7680 

WRITE (66. 30)  CORRD( I TEMP .JTEMP . 1 ) »CQRRD( ITEMP. JTEMP. 2) 

7700 

30 

FORMAT  (IX. 'CENTER  OF  WINDOW  AT  ('.  F9.5' > ' »F?.5» ' ) ' , /) 

7720 

C 

7740 

WRITE  (66.36) 

7760 

36 

FORMAT  (//IX. 'WINDOW  VALUES'./) 

7780 

DO  4  K  =  l.KU 

7800 

DO  5  J  =  l.JU 

7820 

WRITE  (66.6)  (WIND(I»J»K).I=1.IW) 

7340 

6 

FORMAT  ( IX  »<IW>( 14) ) 

7860 

5 

CONTINUE 

7880 

WRITE  (66.8) 

7900 

8 

FORMAT  (/) 

7920 

4 

CONTINUE 

7940 

C 

7960 

WRITE  (66.10)  TIME. XTRANI .YTRANI . MAGI. ROTI .XVELI .YVELI » 

7980 

l 

VMAGI » VRQTI 

8000 

10 

FORMAT  (lX.F5,4.8(3X.E11.4)»/> 

8020 

WRITE  (66.11) 

3040 

11 

FORMAT ( 10X. 'XTRANW' »8X. 'YTRANW' »8X» 'MAGW' » 10X, 'ROTW', 

8060 

i 

10X . ' XVELW ' . 9X . ' YVELW ' . 9X . ' VMAGW ' . 9X . ' VROTW ' ) 

8080 

WRITE  (66.12)  XTRANW. YTRANW. MAGW. ROTW. XVELW. YVELW, VMAGW. VROTU 

8100 

12 

FORMAT  (6X.8(3X.E11,4)./) 

8120 

WRITE  (66.13) 

8140 

13 

FORMAT ( ' 1 ' ) 

3160 

C 

3180 

ELSE 

8200 

WRITE  (66.14)  TIME. XTRANI. YTRANI. MAGI. ROTI. XVELI. TVELI. 

8220 

l 

VMAGI. VROTI 

8240 

14 

FORMAT (1X»F5.4»8<3X»E11.4)) 

3260 

WRITE(66» 15)  XTRANW. YTRANW. MAGW. ROTW. XVELW. YVELW. VMAGW. VROTW 

3280 

15 

FORMAT (11X.8(3X.E11.4>./) 

3300 

ENDIF 

8320 

RETURN 

8340 

END 

8360 

C 

8380 

3400 

c  mm*mmmmmmmt*mmmm******mmt*****m*tm 

3420 

c 

SUBROUTINE  INIT  IS  USED  TO  INITIALIZE  PARAMETERS. 

8440 

c  tmmmmmm*********t***t****ms*****tm****mmm* *** 

9460 

c 

8430 

SUBROUTINE  INIT 

8500 

c 

8520 

c 

*  COMMON  DECLARATIONS  * 

8540 

PARAMETER  (IW=9.JU=9.KW=3.IB=3»JB*3.KB=3.ISIZE=128.JSIZE=128) 

8560 

c 

8580 

COMMON  /WINDOW/  WIND( IW. JW. KW) .CORRD( ISIZE. JSIZE.2) 

8600 

INTEGER  WIND 

3620 

COMMON  /ORIGIN/  10. JO.KO.XO.YO 

3640 

INTEGER  10. JO.KO.XO.YO 
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3660 

COMMON  /FARMS/  ROW, COL, PROW, PCQL, NUMROW, NUMCOL, WTEST 

3630 

INTEGER  ROW, COL, PROW, PCOL, NUMROW, NUMCOL, WTEST 

S700 

COMMON  /COCHGS/  XVELI . YVELI .VROTI .VMAGI .TIME.TSTEP 

8720 

COMMON  /EQU/  ALPHA! 4) . BETA » COEFF ( 9 . 4 ) .VECTOR (9) »SQLN 

8740 

INTEGER  ALPHA. BETA 

3760 

C 

37S0 

C 

THE  PIXEL  VALUE  AT  (COL. ROW)  IS  PUT  INTO  UIND<1. 1,1) 

8800 

WRITE  (6.20) 

8820 

20 

FORMAT  ( IX. 'BEGIN  WINDOW') 

8840 

WRITE  (6.50) 

3860 

50 

FORMAT  (IX. 'ROW  COL') 

8380 

READ  (5.*)  ROW. COL 

3900 

C 

3920 

C 

THE  MOTION  TO  BE  TRACKED 

3940 

WRITE  (6.30) 

3960 

30 

FORMAT  ( IX . '  INPUT.* XVELI  .YVELI .  VROTI/PI .  VMAGI '  > 

3980 

READ  (5.*)  XVELI . YVELI . VRQTI .VMAGI 

9000 

C 

°020 

PROW  =  1 

9040 

PCOL  =  1 

9060 

NUMROW  =  128 

9080 

NUMCOL  =  128 

9100 

c 

9120 

XO  *  COL  +  (IW-l)/2 

9140 

YO  *  ROW  +  (JW-D/2 

9160 

c 

9180 

DO  60  1=1, NUMCOL 

9200 

DO  70  J=l, NUMROW 

9220 

CORRD(I.Jfl)  =  I 

9240 

CQRRD( I, Ji 2)  =  J 

9260 

70 

CONTINUE 

9280 

60 

CONTINUE 

9300 

C 

9320 

TSTEP  =  0.033 

9340 

TIME  =  2*TSTEP 

9360 

PI  =  3,14159 

9380 

VROTI  =  VROTIIPI 

9400 

C 

9420 

10  =  (IW+D/2 

9440 

JO  =  (JW+l)/2 

9460 

KO  =  (KW+D/2 

9480 

C 

9500 

DO  10  1=1,4 

9520 

SOLN(I)  =0.0 

9540 

10 

CONTINUE 

9560 

C 

9530 

FEED(l)  =  1.0 

9600 

FEED<2)  =  1.0 

9620 

FEED<3)  =  1.0 

9640 

FEED<  4)  =  1.0 

9660 

C 

9680 

WRITE  (6,40) 

9700 

40 

FORMAT  (IX, 'INPUT  1  cr.  TO  WRITE  WINDOW  OR  0  cr.  TO  1 

9720 

READ  (5,#)  WTEST 

313 


9740 

9760 

9730 

9800 

9820 

9840 

9860 

9880 

9900 

9920 

9940 

9960 

9980 

10000 

10020 

10040 

10060 

10080 

10100 

10120 

10140 

10160 

10130 

10200 

10220 

10240 

10260 

10230 

10300 

10320 

10340 

10360 

10380 

10400 

10420 

10440 

10460 

10480 

10500 

10520 

10540 

10560 

10580 

10600 

10620 

10640 

10660 

10680 

10700 

10720 

10740 

10760 

10780 

10800 


IF  (UTEST  ,£Q.  0)  THEN 
WRITE  (66,5)  IB 

5  FORMAT  (IX. 'WINDOW  SIZE  =  ',I2,//> 

WRITE  (66. 3) 

3  FORMAT ( IX.' TIME'. 5X. 'XTRANI' »BX» 'YTRANI' .8X.' MAGI' .10X. 'ROTI' 

l  ,10X,'XVELI' »9X,'YVELI' ,9X,'VMAGI' ,9X, 'VROTI' ) 

WRITE  (66.4) 

4  FORMAT ( 15X . ' XTRANW ' . 8X . ' YTRANW ' . 8X . ' MAGH ' . 10X . ' ROTW ' . 

1  lOX.'XVELU' »9X. 'YVELW' »9X» 'VMAGU' .9X. 'VRQTW' , /) 

ENDIF 

C 

RETURN 

END 

C 

c  ******************************************************************** 

C  SUBROUTINE  WINDOW  (1)  USES  THE  INTEGRATION  RESULTS  TO  MOVE  THE 

C  WINDOW  IN  ORDER  TO  TRACK  THE  TARGET  AND  (2)  PERFORMS  THE  MOTION 

C  ON  THE  1ST  IMAGE  USED  TO  GET  THE  WINDOW  AND  WRITES  THE  RESULT 

C  TO  UNIT  =  M  4-  15. 

C  ********************************************************************* 

C 

SUBROUTINE  WINDOW  (M) 

C 

C  *  COMMON  DECLARATIONS  * 

C 

PARAMETER  ( IW=9 . JU=9 »KW=3 . IB=3» JB=3 ,KB=3, ISIZE=128, JSIZE=128> 

C 

COMMON  /WINDOW/  WIND( IW. JW.KW) . CORRD( ISIZE . JSIZE. 2) 

INTEGER  WIND 

COMMON  /EQU/  ALPHA* 4 > , BETA.COEFF ( 9 , 4 > , VECTOR ( 9 ) , S0LN( 4 ) ,FEED( 4 ) 

INTEGER  ALPHA. BETA 

COMMON  /IMAGE/  IMAGE ( ISIZE . JSIZE ) 

INTEGERS  IMAGE 

COMMON  /BUFFER/  BUF ( IB . JB.KB) 

INTEGER  BUF 

COMMON  /FARMS/  ROW, COL. PROW. PCOL. NUMROW. NUMCOL .WTEST 
INTEGER  ROW, COL .PROW, PCOL, NUMROW, NUMCOL, UTEST 
C 

C  t  LOCAL  VARIABLES  * 

LQGICAU1  B IMAGE ( ISIZE  > , BLINE(2*ISIZE,  JSIZE) 

CHARACTER* (ISIZE)  IMAGELINE 
INTEGER*2  NEUIM( ISIZE, JSIZE) 

REAL  ECOSU.ESINW.XU. YW 

EQUIVALENCE  ( BIMAGE, IMAGELINE) ,( BLINE.NEWIM) 

C 

C  AFFINE  TRANSFORMATION  TO  MOVE  WINDOW 

C 

ITEMP  =  COL  -  PCOL  +  (IB+D/2 
JTEMP  =  ROW  -  PROW  +  (JB+15/2 
XU  =  CORRDdTEMP,  JTEMP, 1) 

YW  =  CORRDdTEMP,  JTEMP, 2) 

ECOSW  =  EXP(S0LN(3) )  *  C0S(S0LN(4)) 

E3INW  =  EXP(S0LN(3) )  *  SIN(S0LN(4>) 


j 
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10820 

C 

10340 

DO  10  J=1fNUMR0W 

10360 

DO  20  1=1 fNUMCOL 

10880 

CQRRD< 1 f  J» 1 )=£COSU*(CORRD( I >  J » 1 >-XW)-ESINW*(CORRD( I » J»  2) 

10900 

& 

-YUH-EC0SW*S0LN(1>-ESINW*SQLN<2> 

10920 

CORRD( I f Jf2)=ESINW# (CQRRD( 1 f  Jf  1 ) -XW) +ECOSU*(CQRRD( I » J>2 ) 

10940 

l 

-YW ) 4ES I NW4S0LN < 1 >  +ECOSU*SQLN ( 2 ) 

10960 

CORRD(IfJfl)  =  CORRD( I » J » 1 )  +  XW 

10980 

CORRD( I f Jf 2)  =  YU  +  CORRD( I > Jf2> 

11000 

20 

CONTINUE 

11020 

10 

CONTINUE 

11040 

c 

11060 

c 

OPEN  A  IMAGE  FILE  AND  FILL  PIXEL  VALUES  INTO  THE  ARRAY  IMAGE 

11030 

c 

11100 

DO  30  K=1 f  KU 

11120 

KCOUNT  =  H  H  t  9 

11140 

OPEN(UNIT=KCOUNT i STATUS2 ' OLD' > ACCESS=' DIRECT ' » 

11160 

l 

RECORDTYPE= ' FIXED' »READONLY) 

11180 

c 

11200 

DO  80  1=1 t ISI2E 

11220 

DO  90  J=1 f JSIZE 

11240 

IMAGE ( I >  J )  =  0 

11260 

90 

CONTINUE 

11280 

80 

CONTINUE 

11300 

C 

11320 

DO  40  J=1 .NUMROW 

11340 

READ  (KCOUNT' J)  IMAGELINE 

11360 

DO  50  I=1fISIZE 

11330 

IMAGE(I.J)  =  BIMAGE(I)  .AND.  255 

11400 

50 

CONTINUE 

11420 

40 

CONTINUE 

11440 

C 

11460 

C 

USE  THE  TRANSFORMATION  TO  TRACK?  STORE  THE  RESULT  IN  NEWIM 

11480 

C 

11500 

DO  60  J=1»NUMR0W 

11520 

DO  70  1=1 i NUMCOL 

11540 

IF  (CORRD(IfJfI)  .LT.  1  .OR.  CORRD(IfJfI)  .GT.  isize 

11560 

s 

.OR.  CORRD< If Jf 2)  .LT.  1  .OR.  C0RRD(IfJf2)  .GT.  JSIZE) 

11580 

1 

THEN 

11600 

NEVIIM(IfJ)  =  0 

11620 

ELSE 

11640 

NEWIM(IfJ)  =  IBILIN(C0RRD(IfJf1)fC0RRD(IfJf2)) 

11660 

ENDIF 

11680 

70 

CONTINUE 

11700 

60 

CONTINUE 

11720 

C 

11740 

CLOSE  ( UNIT=KCOUNT ) 

11760 

C 

11780 

C 

FILL  WINDOW  WITH  NEW  PIXEL  VALUES 

11300 

c 

11820 

DO  110  I=1fIU 

11340 

DO  120  J=1fJU 

11860 

ITEMP  =  COL  -  PCOL  +  1 

11380 

JTEMP  =  ROW  -  PROW  f  J 
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11900 

WIND ( I ?  J?  K )  =  NEWIM( ITEMP? JTEMP) 

11920 

1 

20 

CONTINUE 

11940 

110 

CONTINUE 

11960 

C 

11980 

C 

WRITE  TRACKED  IMAGE  INTO  A  FILE 

12000 

C 

12020 

IF  (K  ,£Q.  1)  THEN 

12040 

0P£N(UNIT=KC0UNT415.STATUS=,NEW'»ACCESS='DIRECT' r 

12060 

l  RECORDTYPE=' FIXED' ? RECL=NUMC0L/4?  BLOCKS I ZE=NUMCOL> 

12080 

DO  130  J=1 iNUMRQW 

12100 

DO  140  1  =  1  ?  NUMCOL 

12120 

BIMAGE(I)  =  BLINE(I*2-1?J) 

12140 

140 

CONTINUE 

12160 

WRITE  (KC0UNTH5' JlIHAGELINE 

12180 

130 

CONTINUE 

12200 

CLOSE  (UNIT=KC0UNT+15) 

12220 

ENDIF 

12240 

C 

12260 

30 

CONTINUE 

12280 

RETURN 

12300 

END 

12320 

C 

12340 

C 

ttt%*t**%t%t*tt***t****t*******tt*t*tt*t*t*t*t***t*t*ttt**tt*ttttt ** 

12360 

c 

t*tttt***ttt*ttt*********t*t*tttt*tt*t*t******tttt*t*t*t**t***tt*tt 

12380 

c 

12400 

SUBROUTINE  LINEQ(A?X?B?M?N?CC> 

12420 

INTEGER  CC 

12440 

c 

12460 

c 

SOLVE  AX=B.  T  HOLDS  AN  UPPER  TRIANGULAR  MATRIX  WHILE  S 

12480 

c 

IS  WORKSPACE.  THE  METHOD  FACTORS  A=U*T  WHERE  THE  COLUMNS  OF 

12500 

c 

U  ARE  ORTHOGANAL  AND  T  IS  TRIANGULAR.  THE  RESULTING  SYSTEM 

12520 

c 

T*X=B'  IS  EASILY  SOLVED  BY  BACK  SUBSTITUTION.  ASSUME  M 

12540 

c 

EQUATIONS  AND  N  UNKNOWNS.  (  N  <=  M  <=  9  ) 

12560 

c 

THE  MATRIX  OF  COEFFICIENTS?  A  IS  STORED  IN  THE  FIRST  N  ROWS 

12580 

c 

AND  THE  FIRST  M  COLUMNS  OF  THE  9X9  A  ARRAY.  THE  ROUTINE 

12600 

c 

BRINGS  IN  THE  WHOLE  9X9?  BUT  ONLY  USES  A<1?1>  TO  A(N?M> 

12620 

c 

(RECALL  THAT  FORTRAN  STORES  THE  ARRAY  COLUMN-WISE?  BUT 

12640 

c 

ADRESSES  THE  ELEMENTS  IN  THE  STANDARD  ROW?COLUMN  FORMAT) 

12660 

c 

note:  the  a  array  is  altered  during  execution. 

12630 

c 

12700 

DIMENSION  A(9?9)?T(9?9>?X(N)?B(M) 

12720 

CC  =  1 

12740 

c 

M  MUST  BE  <=  9?  AND  N<=M.  CC  IS  A  COMPLETION  CODE?  IF  THE 

12760 

c 

SUBROUTINE  EXECUTES  PROPERLY  CC  WILL  BE  RESET  TO  0  BEFORE  RETURN 

12730 

DO  40  1=1 ?N 

12800 

IF  (I.EQ.l)  GO  TO  25 

12820 

DO  20  J=1?M 

12840 

3  =  0 

12860 

11=1-1 

12880 

DO  10  K=1  ?  1 1 

12900 

IF  (T <K?K)  .LT.  ,0001)  GO  TO  5000 

12920 

S=S+A(J?K)*T(K?I)/T(K?K) 

12940 

10  CONTINUE 

12960 

A( J? I) =A( J? I ) -S 

316 


i 


12980 

20 

CONTINUE 

13000 

25 

DO  40  K=I.N 

13020 

S=0 

13040 

DO  30  J=1  »M 

13060 

S=S+A( J  . I ) *A( J » K  ) 

13080 

30 

CONTINUE 

13100 

T(I.K)=S 

13120 

40 

CONTINUE 

13140 

DO  60  1=1. N 

13160 

S=0 

13180 

DO  50  J=1 .  M 

13200 

S=S+A< J.I)*B< J) 

13220 

50 

CONTINUE 

13240 

X(I)=S 

13260 

60 

CONTINUE 

13280 

DO  80  1=1. N 

13300 

I 1=N+1 -I 

13320 

IF  (Il.EQ.N)  GO  TO  75 

13340 

12=11+1 

13360 

DO  70  J=I2»N 

13380 

X(I1)=X(I1)-T(I1.J)*X(J) 

13400 

70 

CONTINUE 

13420 

IF  (T(Il.Il).LT. .0001 )  GO  TO  5000 

13440 

75 

X< 11 )=X( 11 )/T( 11 > 11 ) 

13460 

80 

CONTINUE 

13480 

cc=o 

13500 

RETURN 

13520 

5000 

CC=-1 

13540 

C 

A 

COMPLETION  CODE  OF  -1  INDICATES  THAT  THE  SUBROUTINE 

13560 

C 

TRIED  TO  DIVIDE  BY  0. 

13580 

RETURN 

13600 

END 

13620 

C 

13640 

C 

******************************************************************** 

13660 

C 

******************************************************************** 

13680 

c 

13700 

INTEGER  FUNCTION  IBILIN*2(XX.YY) 

13720 

c 

13740 

c 

THIS  FUNCTION  RECEIVES  2  REAL  COORDINATES  (PRODUCED  BY  THE  TRANS¬ 

13760 

c 

FORMATION  IN  THE  CALLING  ROUTINE)  WHICH  ARE  COORDINATES  RELATIVE 

13780 

c 

TO 

THE  OLD  IMAGE.  SINCE  THE  COORDINATES  ARE  REAL  VALUED. 

13300 

c 

THE  POSITION  WILL  NOT  BE  ON  A  PARTICULAR  PIXEL*  BUT  RATHER  AMONG 

13820 

c 

4. 

THIS  FUNCTION  RETURNS  A  BILINEAR  INTERPOLATION  FOR  THE  4 

13840 

c 

SURROUNDING  POINTS. 

13860 

c 

13330 

PARAMETER  (ISIZE=128.JSIZE=128> 

13900 

COMMON  /IMAGE/  IMAGE ( ISIZE .JSIZE ) 

13920 

INTEGER*:  IMAGE 

13940 

REAL  H. V .HTEMP1 . HTEMP2 . VTEMP 

13960 

c 

13980 

MX1=MAX(0. INT (XX) ) 

14000 

MX2=MIN<  ISIZE.  INTCXX)  +■  1 ) 

14020 

MY  1  -MAX ( 0. INT ( YY) ) 

14040 

MY2=MIN( JSIZE. INT(YY)+1) 

i 

< 

i 


I 


'*■  ■  -  -  — * 
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14060  H=XX-rtX 1 

14030  O-YY-rtYl 

14100  C 

14120  HTEMPi =H* IMAGE (MX2  >HY1 ) 4 ( 1 . Q-H) *IMAGE<  NX1 » HY 1 > 

14140  HTEMP2=H*IMAG£(MX2>f1Y2)  +  < 1 . 0-H)*IMAGE(NXl*MY2) 

14160  VTEMP  =V*HTEMP2K1.0-V)*HTEMP1 

14180  IBIIIN*ININT(VTEMP> 

14200  RETURN 

14220  END 

14240  C 

14260  C  tt**t*t**tt*tt**ttt*t*%%***ttt*tttttMtttt*.it*ttt**t**ttt***1<*t*%*t* 
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ADAPTIVE  PATTERN  MATCHING  USING  CONTROL  THEORY  ON  LIE  GROUPS* 


Thomas  3.  Newman  and  Leopold  Zlobec 
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Lubbock,  Texas 


Abstract 

A  method  is  given  for  matching  a  subpattem  of  a  two-dimensional 
image  against  a  stored  prototype,  where  the  latter  is  defined  on  a 
window  whose  position  and  shape  is  determined  by  the  action  of  a  Lie 
group  of  transformations.  The  method  involves  the  construction  of  a 
path  in  the  control  group  along  which  the  matching  error  decreases 
to  a  local  minimum. 


1 .  INTRODUCTION 

A  problem  of  classical  interest  in  pattern 
recognition  is  that  of  determining  the 
presence  or  absence  of  a  particular  sub- 
pattern  or  3ubpattern  class.  In  the  anal¬ 
ysis  of  two-dimensional  imagery  this  can 
take  the  form  of  detection  of  corners  and 
edges  or  the  location  of  a  specific  sil¬ 
houette.  More  particularly,  we  may  be  in¬ 
terested  in  obtaining  an  exact  match  of  a 
specific  portion  of  the  image  to  a  sub¬ 
image,  often  a  prototype,  which  may  appear 
m  an  arbitrary  manner,  varying  in  size, 
location  and  orientation.  This  is  the 
problem  which  is  herein  addressed. 

A  related  question  was  considered  by  Dir- 
ilten  and  Newnan  [3]  where  it  was  shown 


how  two  planar  images  could  be  matched  un¬ 
der  arbitrary  affine  transformation  of  the 
plane,  if  a  match  were  at  all  possible.  In 
addition  to  affine  transformations ,  an  al¬ 
lowance  was  also  made  for  dilation  of  in¬ 
tensity  scale  such  as  that  which  results 
from  under  or  over  exposure  of  film  within 
latitude  limits.  The  results  cited,  how¬ 
ever,  are  of  little  use  in  matching  subpat- 
tams,  since  the  algorithms  are  highly  sen¬ 
sitive  to  the  background  context.  Never¬ 
theless,  the  utility  of  a  group  theoretic 
approach  to  pattern  matching  was  clearly 
demonstrated. 

In  the  following  we  present  a  method  for 
performing  a  local  search  for  an  imbedded 
subpattern  of  a  two-dimensional  image.  The 


'This  research  was  supported  by  the  Army  Research  Office,  Contract 
3AAG29-3C-C-0087  and  by  the  Office  of  Naval  Research,  Contract 
N0014-76-C-1136. 
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method  is  ona  involving  adaptive  control 
of  a  retina  which  seeks  the  desired  sub- 
pattern  by  evolving  along  a  curve  in  the 
space  of  parameters  in  a  direction  which 
assures  improvement  in  the  goodness  of 
fit. 


2 .  BACKGROUND 

Let  G  be  a  Lie  group  of  transformation  on 
an  analytic  manifold  X.  Suppose  G  has  di¬ 
mension  n  while  M  has  dimension  m.  Lat  :< 
and  y  denote  the  coordinates  of  elements  f 
and  g  in  G,  respectively,  in  a  patch  con¬ 
taining  the  identity  element  e  of  G.  Also, 
let  p  denote  coordinates  of  an  element  u 
of  M  in  some  pacch  in  X.  Me  may  then  ex¬ 
press  the  coordinates  z  of  the  product 
h  »  fg  and  the  coordinates  q  of  the  ele¬ 
ment  v  «  gu,  relative  to  suitable  patches, 
by  means  of  analytic  functions 

z  *  J(x,y)  (2.1) 

q  -  K(y,p)  (2.2! 

X  and  J  are  vector-valued,  having  values 
in  n-dimensional  space  Rn  or  Cn  and  a- 
dimensional  space  R3  or  C3.  Hereafter  we 
shall  assume  that  these  underlying  spaces 
are  real,  we  denote  the  ith  component  of 
J  6y  J.  and  the  jth  component  of  X  by  . 

In  order  to  define  the  Lie  algebra  of  G  we 
first  introduce  real-valued  maps  on  G  by 

f>ij(x)  .  (2.3! 

where  i  and  ]  each  range  from  1  to  n.  The 

cross-section  P.^,  which  consists  of  the 
as  i  ranges  from  1  to  n,  and  j  is  fix¬ 
ed,  may  be  thought  of  as  a  vector  field  in 
3  ,  Such  a  vector  field  attaches  to  a 
point  x  the  vector  P.^(x).  As  such,  P#1, 
fora  a  basis  for  the  tangent 
space  at  the  point  x  [1,2) .  The  infinite¬ 
simal  transformations  of  G  may  now  be  de¬ 
fined  by 

n 

:<,  *  '  ?,  < (x)  TT~  -  '2.4! 

i-l  ,xi 


for 


n 


The  differential  operators  so  defined  are 
to  be  considered  as  linear  operators  on  the 
space  of  analytic  functions  on  G,  or,  more 
generally,  on  the  space  of  differentiable 
functions  on  G.  The  Lie  algebra  of  G  is 
simply  the  n-dimensional  vector  space  con¬ 
sisting  of  all  linear  combinations  of  these 
operators,  and  will  be  denoted  iy  KG)  [21. 


The  Lie  algebra  of  G  may  also  be  defined  m 
terms  of  its  actions  on  the  manifold  M. 
.Analogous  to  (2.3!  we  define 


;x 

Q  ;  p!  -  , — — ( y ,  pi  1 
13  3y,  1  &  y-e 

for  1  -  l,2,...,m  and  j  -  1,2,.. 
ally,  as  in  (2.4)  above  we  set 


X' 

3 


i-l 


(2.5) 


?m- 


(2.6) 


The  operators  x£,  x^ ,  ...,  x^  apply  to 
functions  defined  on  M  and  span  a  Lie  alge¬ 
bra  isomorphic  to  L(G). 


The  following  result  from  [41  will  be  used 
later,  and  is  stated  for  reference: 


Theorem  2.1.  Let  f:  M  -  R  be  differenti¬ 
ate  and  define  F:  G  «  M  -  R,  m  terms  of 
coordinates ,  by 

F(x,p>  -  f  (X (x , p) )  .  (2.7! 

Then  for  each  3  -  1,2,..., a  we  have 

X*F  -  X!F.  (2.3) 

j  3 

Let  us  consider  a  curve  t  -  g(t)  in  G  sat¬ 
isfying  g(0)  »  e.  In  terms  of  a  coordinate 
patch  at  a,  git)  may  be  described  by  a 
curve  x(t)  in  Rn  satisfying  x(0)  -0.  we 
shall  consider  the  case  in  which  x(t)  is 
given  as  the  solution  of  an  evolution  equa¬ 
tion  of  the  form 


x  ( t )  -  :  .  (t)P..  (x(tl  )  ,  x  ( 0)  -  3,12.9) 

1  —  1  *  1 

where  *re  cross-sections  of  the 

Array  of  functions  given  by  >2.3) ,  and 
( t „ i t)  are  suitable  control  func¬ 
tions. 


*  ■-■jk 
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Now  let  ?  denote  the  coordinates  of  a 
point  u  in  some  coordinate  patch.  For  a 
differentiate  map  f:  M  *  R  we  may  define 
H :  a  <  M  -  R  by  setting 

H(t,p)  -  f (gCt) u) .  (2.10) 

we  recognize  that  H(t,p)  »  F(x(t) ,p)  where 
F  is  the  extension  of  f  to  G  <  M  as  in 
Theorem  2.1  above.  From  the  point  of  view 
of  application,  if  we  regard  f:  M  -*  R  as 
an  image,  then  H(t,p)  represents  the  mov¬ 
ing  image  obtained  by  translation  due  to 
the  curve  g(t).  Also  from  (41,  we  have 

Theorem  2.2.  In  the  context  above , 

3  V  - 

ytr  ™  1,  (t)X!H.  (2.11) 

Jt  iml  i  i 


3.  THE  CONTROL  MODEL 

3y  an  image  we  mean  a  map  f:  M  -  R,  where 
the  value  f(p)  at  a  point  p  e  M  represents 
the  gray  value  at  the  picture  element  at 
p.  In  practice,  values  are  observed  on  a 
subset  W  M,  which  we  regard  as  a  window 
which  may  be  translated  by  the  action  of  G 
on  M.  Thua,  upon  translation  by  an  ele¬ 
ment  x  e  G,  the  value  observed  at  p  6  W  i3 
given  by  F(x,pi  «  t (Xlx. p) ) .  as  in  (2.7) 
above. 


We  consider  a  given  prototype  sub- image  v 
defined  on  the  window  W,  V:  w  -  R.  The 
problem  then  is  to  determine  x  S  G  such 
that  F(x,p)  ■  V(p)  for  all  p  e  w,  or  deter¬ 
mine  that  no  such  x  exists.  As  a  matter 
of  practice,  we  seek  x  6  G  which  minimizes 
the  objective  function 


?(X)  -  J  f  (F{x,p)  -  V(p))2dp,  (3.1) 

W 

where  dp  represents  a  volume  element  and 
the  integral  is  over  the  window  w,  which 
is  assumed  to  be  of  bounded  volume. 

In  general,  for  any  two  functions  f^,f, : 
w  -  M  we  define 


<w 


w 


and 


1  :f. 


<!v!i 


.1/2 


Thus,  • (x)  3  F  -  V  ~ / 2 ,  where  x  is  re¬ 
garded  as  a  parameter. 

The  following  is  a  well-known  property  of 
the  Lie  group  G  (2): 

lemma  1 .  In  order  that  the  differential 

d  (x)  *  0  at  a  point  x  S  G,  it  is  necessary 

and  sufficient  that  each  X;?(x)  •  3  where 

X, ,X,,...,X  are  the  generators  of  L(G) 
i  4-  n 

given  by  (2.4)  . 

3y  direct  calculation,  we  obtain  X,?(x)  » 

/(F!x,p)  -  V(p) ) X; F (x,p) dp.  In  practice, 

W 

this  expression  is  difficult  to  compute 
numerically,  due  to  the  presence  of  the 
term  X^F,  which  cannot  he  computed  directly 
from  observed  data.  However,  by  Theorem 
(2.1)  we  have  x^F  »  xj^F,  and  the  latter  can 
be  calculated  from  a  single  value  of  x. 

Suppose  now  that  a  curve  in  G  is  given  by 
coordinates  x(t)  obtained  as  a  solution  of 
Equation  (2.9).  we  seek  to  find  i ( t)  • 
(\1(t)  , . . .  .'jjte) )  so  that  silt)  »  'Mx(t)) 
decreases  to  a  minimum  value.  Defining 
H ( t , p)  •  F (x ( t) ,p)  we  obtain, 

L(t)  *  /  ( H  ( t ,  p)  -  V(p))i£L  (t,p)dp  (3.2) 

which,  by  application  of  Theorem  (2.2),  be¬ 
comes 


■Mt)  •  I  \.  (t)  ;  (H(z,p)  -  V(p))x;H(t.p)dp 
i-1  1  /.  1 

n  <3-3> 

-  i  > .  (t) <h  -  v,x:h> 

i«l  “  1 


Upon  observing  that  <H  -  V,.\)H>  » 

<F  -  V,XjF>  *  x^f  at  x  »  x(t),  we  deduce: 

Theorem  3.1.  If  i^lt)  is  chosen  so  that 
sgnAi(t)  *  -  sgn  <H  -  v,x^H>,  we  have 

p(tl  i  0  for  all  c,  with  equality  at  t  » 
if  and  only  if  i’t  »  0  at  x  «  x(tg)  • 

Among  the  claas  of  bounded  controls, 

' -  (t)  i  1,  we  see  that  the  rats  of  de¬ 
crease  of  (t)  is  maximized  by  the  choice 
• L ( t )  *  -  sgn  <H  -  7,X!H',  (3.4) 


i 
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for  i  ■  1,2 . n.  Of  course,  ocher  stra¬ 

tegies  can  be  formulated,  including  steep¬ 
est  descent,  and  some  methods  using  un- 
oounded  controls.  3y  proceeding  along 
trajectories  defined  by  the  solution  of 
,2.9!  with  '(t)  given  by  (3.4),  we  ap¬ 
proach  a  critical  point  of  *  ;i.s.  d  *  0). 
Since  maxima  and  saddle  points  are  un¬ 
stable  . nder  perturbation,  in  practice 
this  extrema  point  will  always  be  a  mini- 


4 .  SIMULATION  RESULTS 

The  results  discussed  in  the  previous  sec¬ 
tion  have  been  implemented  by  a  discrete 
algorithm  and  tested  on  simulated  data  ( S  ]. 
A  digitized  two-dimensional  image  was 
first  generated  in  the  form  of  a  large 
two-dimensional  array,  and  the  prototype 
was  generated  it  a  20  ■  20  window  array. 

The  image  space  was  assumed  to  be  subject 
to  translation,  magnification  and  rotation, 
giving  rise  to  a  four  parameter  Lie  group 
of  transformations  in  the  plane,  R2. 

A  number  of  cases  were  considered,  includ¬ 
ing  sore  involving  multiple  (false)  tar¬ 
gets  and  others  in  which  the  prototype  was 
absent  from  the  image  being  searched.  In 
some  cases  the  image  was  contaminated  by 
3%  random  noise,  rn  all  cases  the  search 
was  started  with  overlap  between  the  pro¬ 
totype  target  and  the  image  target. 

The  differential  equation  (2.9)  was  solved 
by  means  of  a  Runge-Kutta  fourth  order 
method,  with  a  dynamic  step  size,  which 
was  increased  as  necessary  to  accelerate 
convergence  and  decreased  as  necessary  to 
maintain  stability.  Integration  was  re¬ 
placed  by  summation,  although  we  conjec¬ 
ture  that  convergence  could  have  been  ac¬ 
celerated  by  the  use  of  a  trapezoid  rule. 

generally,  search  times  ranged  from  30  to 
30  steps,  with  the  longer  search  times 
prevailing  for  the  more  difficult  cases. 


In  ill  cases,  the  final  results  ware  quite 
reasonable ,  even  m  those  cases  where  the 
prototype  was  absent.  In  the  latter  cases, 
the  search  terminated  with  a  "best"  match, 
with  a  commensurately  large  final  error. 

As  an  example.  Figure  1  shows  that  starting 
position  for  a  noisy  image  containing  two 
objects.  The  prototype  is  indicated  by  the 
central  silhouette,  while  the  true  target 
is  shifted  upward,  slightly  to  the  right 
and  is  reduced  m  size.  A  false  target 
overlaps  the  lower  right  comer  of  the  pro¬ 
totype  . 


Fig.  1.  Initial  Window  Position. 

The  termination  conditions  are  shown  in 
Figure  2,  where  the  true  target  was  located 
after  49  steps.  All  parameters  were  cor¬ 
rect  with  the  exception  of  magnification, 
which  was  about  51  too  large.  Smaller  val¬ 
ues  of  magnification,  however,  increase  the 
error  due  to  the  presence  of  the  false  ob¬ 
ject,  which  is  barely  touching  the  bottom 
edge  of  the  window  m  Figure  2. 
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Fig.  2.  Terminal  Window  Position. 
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